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ABSTRACT 


Analytical modelling of heterogeneous and, especially, fractured reservoir systems 
constitutes an important aspect of petroleum reservoir engineering and formation evaluation 
and environmental engineering studies. Fractal geometry has been used in some recent 
studies to help define hydrologic conceptual models for a certain class of fractured systems. 
However, transient pressure behaviour during the flow of non-Newtonian power-law fluids 
in homogeneous systems has often been misinterpreted as flow in fractures when using 
Newtonian-based theory. Thus, use of the existing fractal pressure transient models may not 
be appropriate when analyzing pressure transient data during non-Newtonian fluid flow 


through a network of fractures. 


The principal objective of the present study is to examine the basic characteristics of 
transient flow of non-Newtonian power-law fluids in a fractal network of fractures. The 
theoretical basis for analyzing transient pressure and rate data under such situations are 
presented. Analytical and finite-difference solutions of the partial differential equations 
describing the single-phase flow of a slightly compressible, power-law fluid in an infinitely 


large fractal reservoir are obtained. 


Transient pressure behaviour for finite reservoir cases is studied. Analytical solutions 
are also presented for the case of a two-zone composite reservoir, with the inner zone being 
completely fractal and the outer zone homogeneous. Finally, an analytical model is presented 
for Newtonian fluid flow in a "double-porosity" situation where the fracture network is 
characterized by a particular type of fractal geometry and the matrix is homogeneous. The 
sensitivity of the system response to the various relevant parameters is also examined in 


detail. 
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a percolation system 

dimensionless time defined by Equation (10-8) 

dimensionless time defined by Equation (4-35) 

dimensionless skin factor 

superficial velocity in the radial direction [m/s] 

dimensionless parameter defined by Equation (5-21) 
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dimensionless parameter defined by Equation (8-30) 
variable defined by Equation (6-1) 

parameter defined by Equation (5-17); also group defined by 
Equation (6-18) 

parameter defined by Equation (8-26) 

parameter defined by Equation (5-18) 

parameter defined by Equation (8-27) 

parameter defined by Equation (5-19) 

parameter defined by Equation (8-28) 

gamma function 

incomplete gamma function 

parameter defined by Equation (8-38) 

pressure drop [Pa] 

pressure drop at t = 0 [Pa] 

skin pressure drop given by Equation (5-61) [Pa] 
conductivity scaling exponents of a percolation system for Newtonian 
flows 

conductivity index (fractal exponent) 

matrix/fracture interaction parameter 

viscosity [Pa.s] 

apparent viscosity [Pa.s] 

"effective viscosity" defined by Eqaution (4-16) [Pa.s®.m!-1] 
density [Kg.m-3]; dimensionless variable defined by Equation (5-13) 
dimensionless variable defined by Equation (8-24) 
fracture/matrix interaction index 

dimensionless group defined by Equation (8-7) 
dimensionless group defined by Equation (8-14) 

porosity 

fracture porosity 

matrix porosity 

dimensionless storage parameter 

dimensionless time defined by Equation (9-12) 

group defined by Equation (9-19) 
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Chapter I 


INTRODUCTION 


Research in modelling of pressure transient behaviour of fractured hydrocarbon 
reservoirs has advanced considerably over the past few decades. In 1963, the now classical 
study of Warren and Root's two-scale model of naturally fractured reservoirs was published. 
Since then, most conceptual models of flow in fissured systems have generally centred on 
various forms and modifications of the original double-porosity model envisaged by 
Warren and Root (1963). These conventionally accepted models for interpretation of 
transient pressure data from naturally fractured reservoirs usually involve simplistic 
assumptions regarding the geometry and transport behaviour of fracture networks. More 
importantly, pressure transient responses predicted by these models are sometimes not 


observed in actual well tests in naturally fractured reservoirs. 


In a few recent studies (Chang and Yortsos, 1990; Beier, 1990a; Beier, 1990b; Acufia 
et al., 1992a; Acufia et al., 1992b), the concept of fractal geometry has been made use of in 
order to develop new models for interpretation of pressure transient tests in fractured media. 
In these conceptual models, fractal properties have been attributed to networks of fractures in 
fissured rocks so that the hydraulic response of the fracture system can be analyzed. The 
application of fractal geometry to the analysis of pressure transient behaviour was the result 
of an extension of the work of physicists who studied diffusion in disordered media and 
fractal objects (Orbach, 1986; O'Shaughnessy and Procaccia, 1985). Since pressure for fluid 
flow in porous media satisfies a diffusion-type equation, scaling principles for diffusion in 
fractal objects were applied by analogy to fluid flow through porous media (Beier, 1990; 


Acufia et al., 1992a). 
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The general theoretical formalism of the application of fractal geometry to pressure 
transient analysis was presented by Chang and Yortsos (1990). Their approach was applied 
to analyze real well-test data in a few subsequent studies (Beier, 1990a; Acufia er al., 1992a; 
Acuna et al., 1992b). However, the analysis presented in these studies is restricted as it 
applies specifically to the transient flow of a Newtonian fluid in an infinite flow system. The 
transient injection (or production) behaviour of non-Newtonian fluids is an important and a 
frequently encountered phenomenon for the petroleum reservoir engineering community, 
especially during the flow of fluids such as polymer solutions, emulsions, fracturing fluids 
and oil-sand/fines mixtures through porous media. Without the aid of a proper tool for 
modelling transient flow behaviour of non-Newtonian fluids, analysis of pressure transient 
data during the flow of such fluids would be difficult. Moreover, it is to be noted that the 
results of Chang and Yortsos (1990) apply only for infinite reservoirs and that they also 
assume that the reservoir exhibits fractal characteristics over all length scales. For reservoirs 
exhibiting fractal characteristics only over a finite region around the wellbore, alternative 


approaches must be sought. 


The purpose of this investigation is to address some of the concerns mentioned above. 
Broadly, this work deals with the development of an approach for the interpretation of 
pressure transient tests during the flow of a non-Newtonian power-law fluid in a totally 
fractal reservoir. Specifically, this study examines the pressure transient response of a 
single-well system in a manner consistent with the expectation that the fracture network 


dominates the flow behaviour in the naturally fractured reservoir. 
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Chapter II 


REVIEW OF LITERATURE 


The basic aim of this study is to develop a pressure transient model for the flow of a 
single-phase non-Newtonian power-law fluid through a fractal reservoir. Thus, the topic to 
be discussed in this research can be described broadly as the incorporation of the rheology of 
power-law fluids into a pressure transient model for a fractal reservoir. Accordingly, the 
literature review has been divided into two sections: pressure transient modelling of power- 


law fluid flow, and fractal pressure transient analysis. 


2.1 Literature Survey of Pressure Transient Models for Power-Law Fluid Flow 


The pressure transient modelling and analysis effort in the petroleum industry has 
generally concentrated on the behaviour of Newtonian fluids in homogeneous and 
heterogeneous reservoirs. Considerable numbers of exact and approximate models and type- 
curves exist for the analysis and interpretation of well-test data for various kinds of reservoir 
systems and wellbore conditions and for different tests. However, not much work or data is 
available in the petroleum engineering literature on pressure transient modelling of non- 
Newtonian fluid flow in porous media. In fact, some of the most popular textbooks on well 
testing (Earlougher, 1977; Streltsova, 1988; Sabet, 1991) do not consider non-Newtonian 
pressure transient analysis. Most polymer solutions, emulsions (Olarewaju, 1992), aqueous 
foams (Ikoku, 1978), solid particles suspended in Newtonian liquids (Barnes, 1989) such as 
heavy oil mixed with fines (Poon and Kisman, 1992), hydraulic fracturing fluids (Torok and 
Advani, 1987), and so forth, are non-Newtonian. Approximation of the pressure transient 
behaviour of these fluids by Newtonian fluid flow models, as is commonly done 


(Olarewaju, 1992), may result in significant errors in the ensuing analysis. Thus, from a 
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petroleum reservoir engineering point of view, there is a need for a thorough inspection of 


the transient flow of non-Newtonian fluids in porous media. 


A number of articles on the flow behaviour of non-Newtonian fluids through porous 
media have been published in the chemical engineering, rheology and petroleum engineering 
literatures. In 1969, Savins presented a review of contemporary literature on the rheological 
behaviour of non-Newtonian fluids flowing through porous media. He also discussed the 


relevance of non-Newtonian flow through porous media to various technology areas. 


McKinley et al. (1966) studied the linear flow of a polymer solution in porous media. 
They proposed a modification of the viscosity to be used in the Newtonian Darcy equation 
to model linear flow of a non-Newtonian fluid. They suggested that the viscosity should not 
only be expressed as a function of the rheological properties of the fluid, but also as a 
function of the characteristics of the porous medium and the pressure gradient. Gogarty 
(1967) extended the work of McKinley er al. (1966) and suggested that the effective 


viscosity should also be related to the average shear rate in the core. 


One of the earliest studies of transient pressure behaviour during the flow of a non- 
Newtonian power-law fluid in a porous medium was presented by van Poollen and Jargon 
(1969). They proposed a finite-difference model of a radial system to predict the transient 
flow behaviour of power-law fluids. Their model made use of Newtonian fluid flow 
equations, with the non-Newtonian effects being incorporated by means of the viscosity 
varying as a function of the radial location. A finite-difference model was also presented by 
Bondor et al. (1972) to simulate polymer flooding. However, transient flow was not 


considered in this study. 


Ikoku and Ramey (1979) and Odeh and Yang (1979) presented analytical models for 


transient flow of non-Newtonian power-law fluids in infinite reservoirs. In these studies, 
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however, the effect of wellbore storage was not considered. In a subsequent study, Ikoku 
and Ramey (1980) extended their previous model (1979) to flow in finite circular reservoirs. 
They also used a wellbore storage simulator to study the effects of skin and wellbore storage 


during the transient flow of power-law fluids in infinite and finite systems. 


Pascal and Pascal (1985) studied the steady and unsteady flow of power-law fluids 
through porous media. In order to deal with the nonlinearity of the equations describing the 
transient flow of power-law fluids, they used generalized Boltzman transformations for 
linear and radial flows. They also discussed the limitations associated with the approximate 
analytical solutions of the transient power-law flow problem derived by Ikoku and Ramey 


a 


In a recent study, Vongvuthipornchai and Raghavan (1987) examined numerically the 
pressure falloff behaviour in fractured wells after the injection of a non-Newtonian power- 
law fluid. The problem was studied previously by Murtha and Ertekin (1983); Murtha and 


Ertekin, however, did not present a method to analyze the pressure falloff data. 


Recently, Olarewaju (1992) presented a study to demonstrate the difference in transient 
behaviour between a Newtonian and a non-Newtonian power-law type fluid in 
homogeneous and double-porosity reservoirs. He presented his mathematical solutions to 


develop pressure and pressure-derivative type curves for such reservoir systems. 


2.2 Literature Survey of Fractal Pressure Transient Analysis 


The traditional double porosity model (Warren and Root, 1963) has been considered 
as the standard tool for analysis of pressure transient tests in naturally fractured reservoirs. 
Essentially, this model assumes that a fractured formation consists of two coexistent 


components characterized by two distinct permeability and porosity scales. The fracture 
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system is assumed to be homogeneous (characterized by a Euclidean dimension, dg e.g., dg 
= 1 for a single fracture) and is embedded within the matrix which is also Euclidean (e.g., of 
embedding dimension d = 2 for a cylindrical symmetry flow system). An extension of such 
models was proposed by Abdassah and Ershaghi (1986) who presented a triple-porosity 
model for analysis of pressure transient behaviour of fractured media. This model assumes 
that fractures have homogeneous properties throughout and interact with two groups of 
matrix blocks having different porosities and permeabilities. The approach of Abdassah and 
Ershaghi (1986), however, is less applicable to systems where the fracture network is not of 
Euclidean geometry (i.e., the fracture system is non space-filling) (Chang and Yortsos, 
1990; Acufia et al., 1992a). In a few recent theoretical studies, several authors (Chang and 
Yortsos, 1990; Beier, 1990a; Beier, 1990b; Acufia and Yortsos, 1991; Acufia et al., 1992a; 
Acufia et al., 1992b) have proposed the concept of fractal geometry as an alternative to the 


classical analysis of pressure transient data for various naturally fractured reservoirs. 


The term "fractal" has been used by researchers to characterize geometric objects 
whose properties show a power-law type spatial dependence (Chang and Yortsos, 1990; 
Beier, 1990a; Beier, 1990b; Acufia et al., 1992a; Acufia et al., 1992b). Acceptance of a 
fractal model within the drainage area of a well implies porosity and permeability 
distributions that exhibit power-law dependence on the distance, r, from the well. Chang and 
Yortsos (1990) presented a model for pressure transient analysis of fractal reservoirs where 
the average porosity and permeability of the fracture network over a region of scale r vary 
with r in a power-law manner. They examined the unsteady-state flow of a slightly 
compressible fluid in a fractal fracture network embedded in a Euclidean matrix. An 
appropriate modification of the diffusivity equation for such a case was also undertaken. 
Solutions of this formulation were obtained for the case where the matrix does not 
participate in the flow process and also for the case where the matrix exchanges fluid with 


the fracture network. 


1.8) viet i , ean . 
rs sei ays ot pati gh halle at Se mds | 
(age panne Cision i ta ats ae 
sestitniur Eebaoet alt tea tieee sai Se stort iesseoxa mune i iets ua mia it 
co Sage wet juve iglssind onstage saesatiisin, # anemquenad wed s) 


5 a 
ee ees) Ls fered | wei whe es wie Naas SB Li ens teh: iget, ww ABR 


U _ Mae #4. # i e 
aAatd S SP ie v¢ : pial ae 7 
wy on ty teh ee ee a sealer ea niet oil Ohaaainy ied et aware Oe > ia 
aS Snr ghey (oni agt os Ut aeetEy? sissd] oS O16! CUS 
tury eppegt Dy a Age sane Hee TOD ate ast erie mne eer Wert eal fet ef a ta ae 
f : Lae | 
tyes ata CIRC? awl hte ieek oe ith e148 RUM is IA 
, 4 ' hae Ap ae FR Seine nat teh Te Mgte Nias mie pee Le ie ? 7a! we Loat a ro OF 
, A — 
iy ae Lehirest wuetiie a Ae WT) aay ee eee He, \o % oo 
idee Sone dares 2h i ware ed Dae @eed ret “he ar)" eom-gag 
oun 
(CA) acure t Sne oie raorgit a Wii o> athens <é woe ee 


ot Saas Jde a A ade CYT OS weeny a ORY | ist MARY . 

; ; 7 Wass = 
igen bawepy Maye iy) Som Ocune <3 anlw 2 z 
toe rear. Vas: cohen septuiesithy al cata wthebenaiggts Wal tyre sietlic ey ORE sayTetl 


: v ae 


ube cee en abayand eerie ty a yor aoe e hwaserny 7, 


¥ rel, eee ‘ Per et 7) 


a 6 


thy Rene hy ieee etal Potion, sade Wo -fiicianereny Oe ‘erica Ce = 


. or 4 
eo 


Aaa ee well ancient Paulas vedi sana walaoeng st 


7 


Hh opvtuab eractihit aim ee cepraesiap nieith ‘SI 
"raider may ceglareie weet @ Mattei, hau hese em ‘hits 


* MOTE 


= a ng ml on vy 
4 7) 
ee 


ae “a : 


Beier (1990a) presented an extension of the work by Chang and Yortsos (1990). 
Beier's model applies specifically to Newtonian flow in a cylindrical symmetry reservoir (d 
= 2) containing a fractal permeable network. Beier also chose to write the pressure transient 
equations for a fractal reservoir in a form that required available estimates of "near wellbore 
porosity and permeability". With a proper choice of dimensionless variables, it can be 
shown that the dimensionless pressure transient equation derived by Beier (1990a) is quite 
similar to one of the equations solved by Chang and Yortsos (1990). The main difference in 
the formulation of these equations lies in the way the dimensionless variables were defined 
in these two studies. Beier (1990a) also exhibited some field data from the Grayburg and 
San Andres formations in southeastern New Mexico that apparently do not match solutions 
from conventional models. Instead, Beier's fractal reservoir model was able to provide a 
quantitative analysis of the field data. Similar cases have also been observed from some 
west Texas reservoirs! where the use of a fractal reservoir model provided a much superior 
match to the field pressure data than that provided by any other existing pressure transient 
models. In a subsequent study, Beier (1990b) presented a model to analyze the pressure 
response of a well with a vertical fracture in an infinite fractal reservoir. He showed that 
during the early linear flow period in such a system, the slopes of the log-log plots of 
pressure and pressure-derivative curves are greater than one-half. Because of negligible 
wellbore storage effects in these tests, he argued, homogeneous reservoir models could not 


explain such "anomalous" transient pressure behaviour. 


In 1988, Barker presented his "generalized radial flow model" which has since then 
been of considerable interest and practical use to researchers in hydrogeology and 
environmental engineering. Barker developed his model for transient flow during hydraulic 


tests in fractured media having homogeneous and isotropic properties. In Barker's approach 


TR. A. Beier: private communication, 1992. 
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(1988), the generalized flow equations are developed by assuming that the areas (open to 
flow) of surfaces perpendicular to flow vary in a power-law fashion with distance from the 
centre. Theoretically, this is equivalent to assuming that conductivity and storativity have 
power-law dependence on the radius with the same exponent (Acufia et al., 1992a), so that 
the diffusivity has no spatial dependence. In a subsequent study, Doe (1991) extended 
Barker's approach (1988) by considering unusual shapes of drainage volumes which can 
give rise to power-law variation of flow area with radius. Doe also presented an application 
of Barker's fractional (or generalized) dimension theory to constant pressure tests. Barker's 
fractional dimension model has actually been used to match successfully various field 
hydraulic test data, which could not have been matched by any conventional pressure 
transient models”. Barker's model has also been used in analyzing pressure-transient 


behaviour of nuclear waste repositories in fractured granite (Long et al., 1990). 


The theoretical model of Chang and Yortsos (1990), for the case of a line-source well 
producing from an infinite medium, was used by Acufia et al. (1992a) to interpret the fractal 
characteristics of a naturally fractured geothermal field. In a more recent study, Acuna et al. 
(1992b) reviewed the theoretical background of fractal analysis. They also demonstrated the 
application of various diagnostic techniques for fractal pressure transient analysis as 
developed by Chang and Yortsos (1990). They presented a discussion of various large-time 
behaviour of the pressure and pressure-derivative responses for a fractal system. The 
authors also commented on the implications of a transition in flow (fractal) dimensionality 
during a well test in a fractal reservoir. For example, the early-time response of a well test 
(in a three-dimensional space) may be characterized by a flow dimensionality indicative of 
fluid flow not only in the areal plane but also in the vertical direction (caused by, say, a 


partially-penetrating well in a thick formation); at larger times, however, the pressure 


2 T. W. Doe: private communication, 1992. 
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response may be more likely to be characteristic of cylindrical flow (due to a 
small/negligible vertical flow component), thereby indicating a different dimensionality. This 
transition may impart a double porosity-like "V" shape in the pressure-derivative behaviour, 


causing an incorrect interpretation of the system response. 


Chapter III 


OBJECTIVES AND APPLICATIONS 


3.1 Purpose and Scope of Present Study 


Pressure transient analysis is one of the primary tools used by the petroleum reservoir 
engineering community to characterize the conductive and storage properties of a reservoir. 
The techniques used to analyze the pressure data collected during a period of 
production/injection are based on solutions to partial differential equations describing the 
flow of fluids through porous media. A major assumption incorporated into these solutions 
is that the conductive and storage properties are uniform in the unit of interest, i.e., the 
properties are invariant in space. Moreover, the majority of the existing pressure transient 
models assume Newtonian fluid flow in reservoirs. Although these assumptions have been 
shown to be viable in many situations, the geological complexity of some units and/or the 
non-Newtonian characteristics of certain fluids may make these assumptions of dubious 
validity. More importantly, it has been observed in some previous studies (Odeh and Yang, 
1979; Olarewaju, 1992) that pressure data collected during the flow of non-Newtonian 
power-law fluids may show anomalies when analyzed using methods derived for 
Newtonian fluids. Such anomalies usually have appearances similar to those for flow in 
fractures. And as will also be demonstrated in the present study, pressure behaviour during 
the flow of a power-law fluid in a homogeneous system may appear quite similar to that for 
Newtonian fluid flow in a fractal reservoir. Thus, use of existing pressure transient models 
may not be appropriate for the analysis of pressure data for non-Newtonian power-law fluid 
flow in fractal reservoirs. The purpose of this study is to examine methods for the analysis 
of non-Newtonian fluid flow data in one class of geological settings that is not amenable to 


the conventional approach. 
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In the present study, an investigation of the radial flow of non-Newtonian, power-law, 
slightly compressible fluids in heterogeneous cylindrical-symmetry reservoirs with a fractal 
structure is undertaken. The main objectives of this work are to derive a new partial 
differential equation describing power-law flows in fractal reservoirs and to obtain analytical 
solutions of this equation for a finite-wellbore case in a single-well situation. Examination of 
the system response is conducted under the assumption that the fracture network dominates 
the flow behaviour. Steady and unsteady flow situations in an infinite system are considered. 
Transient pressure behaviour in finite-sized systems are also dealt with. Analysis of the 
system response when the matrix participates in the flow is performed for a simple limiting 
case. Finally, some comments are made on the application of the model to a flow system 
comprising a two-zone composite reservoir, with the inner zone being totally fractal and the 


outer zone homogeneous. 


3.2 Potential Applications 


This study, in effect, extends the modelling effort of non-Newtonian fluid flow to 
behaviour in fractal reservoirs which so far has been confined to Newtonian flow. The 
results of this study should prove to be useful to researchers in the areas of petroleum 
reservoir engineering, hydrogeology and environmental engineering where fracture flow 
modelling is a subject of considerable theoretical and practical interest. The present model 
may be used to interpret transient flow of polymer solutions, emulsions and aqueous foams 
through areally heterogeneous porous media with a fractal structure. The model can be 
applied also to heavy oil reservoirs, like Lloydminster in Alberta, where the reservoir fluid 
of heavy oil mixed with sand has been modelled as a non-Newtonian power-law fluid (Poon 
and Kisman, 1992). The heterogeneous nature of this reservoir may render it to be a 


probable candidate for representation by fractal geometry. 
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Chapter IV 


MATHEMATICAL MODEL FOR POWER-LAW FLOW 
THROUGH FRACTAL POROUS MEDIA 


In this chapter, the rheological properties of non-Newtonian fluids and physical 
characteristics of fractal permeable networks will be considered first to set the conditions of 
the present study. Various terms related to the rheology of non-Newtonian fluids and fractal 
reservoirs will be defined. Subsequently, we will discuss the derivation of a nonlinear partial 
differential equation describing single-phase flow of a non-Newtonian power-law fluid in a 
reservoir showing fractal structure. Finally, we will show the use of a "linearizing 
approximation" that reduces the exact nonlinear equation to a linear one for which closed- 


form analytical solutions may be obtained. 
4.1 Flow of Non-Newtonian Power-Law Fluids Through Porous Media 
4.1.1 Definitions 


A knowledge of the viscosity of various fluids is essential for a study of the flow 
behaviour of these fluids through porous media. The viscosity of a fluid, a measure of its 
resistance to flow, is determined by the transport of momentum in a direction perpendicular 


to the direction of flow. In the viscous region for the flow of a Newtonian fluid, viscosity is 


defined by 
Sih de icy 41) 
dy 
Or 
pee (4-2) 
y 


12 


wey er case eee 
7.4 x a. 7 7 - 
‘Wee Hwa aa 
‘sc 


las wenn pst shia a Tes ante: i tig “4 pao a Je wade 
aia SeOre sTeiwrebaehing at shorn sh siecle 
4 Oye ath i tasarer Op the eojabcndnth Sy feeesiee it spac cere on 


bs pong tasniticvtsh @ 10 Odi an ag Mt oily iter cw (almaapence dante ed jie ie ‘ 


Hy Mol gree citar fy wet; Me inal fh thi chgte aieiben wore 


ch 
runt” BAG ae 5 uP vot Wn Pe Tan iy rious ident iio tt 7 
hos ‘ ite | ‘oie vl Mie o?) 4 oer OMe ad ‘en ae Saeyesirk glee ats gst vadir "s il 7 


| iia dam cocstloe Lenk cha 


ribabs Kay tat Ndee “et aly! rises Leeann ane 


@ 1. ; 


i 


sists ats exp A allan 3 (nartee aa) abyiallt hate os ‘nay Yo oS 
a tosssbtbalebe ¢ ey a 0 pinay ait afar vunnen Pago . abdul nad 
yates ene aehoel hg tg dr 
i leacait 2h cians hapten denne wt ad 
et stn 


where j/ is the Newtonian viscosity, Tt is the shear Stress in the shear plane parallel to the 
flow direction and y is the shear rate perpendicular to the plane of shear. Fluids that obey 
Equation (4-1) are termed "Newtonian fluids". The connotations of the term "Newtonian" 
are (Harris, 1977): 

1) the stress acting on an element of material is proportional to the shear rate, with the 
constant of proportionality being called the viscosity; 

2) the viscosity is time-independent and is also independent of shear rate and time- 


derivatives or integrals of shear rate to any order. 


The exceptions to Newton's viscous law are not rare. There are fluids like polymer 
solutions, colloids, foams, solid-liquid suspensions, etc., that do not exhibit Newtonian 
behaviour. These fluids, called "non-Newtonian fluids”, do not show a direct proportionality 
between the shear rate and the applied shear stress. The viscosity of these fluids changes 
significantly with shear rate under isothermal conditions. It is often convenient to define an 


"apparent viscosity" for non-Newtonian fluids as 


T 
Happ = ¥ (4-3) 


which is a function of shear rate. The rheological properties of time-independent non- 
Newtonian fluids depend only upon the magnitude of the shear stress and not upon the 
duration of the stress. Shear-thinning (pseudoplastic), shear-thickening (dilatant) and 
Bingham plastic fluids fall into this category. With increasing shear rate, pseudoplastic and 
dilatant fluids exhibit decreasing and increasing apparent viscosities, respectively. Bingham 
plastic fluids are those for which a finite shearing stress is required to initiate motion and for 
which there exists a linear relationship (beyond the yield stress) between shear stress and 


shear rate. 
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4.1.2 The Ostwald-de Waele Power-Law Model 


For the purposes of modelling and calculation, the most commonly implemented 
rheological model used to describe shear-thinning or -thickening behaviour is the Ostwald- 
de Waele power-law model. This two-parameter function which has been useful in fitting 
rheological data for a large variety of shear-thickening and shear-thinning flows is often 
represented by the following (Barnes, 1989; Harris, 1977; Schowalter, 1978; Ikoku and 
Ramey, 1979; Torok and Advani, 1987; Poon and Kisman, 1992; Olarewaju, 1992): 


GS MEGe (4-4) 


where H is the consistency and n is the dimensionless flow behaviour index. Equation (4-4) 
"is probably the ... most widely used equation in all of rheology” (Schowalter, 1978). The 
appeal of the power-law is evident. When n = 1, Equation (4-4) reduces to the description of 
a Newtonian fluid with viscosity H. For 0 <n < 1, Equation (4-4) describes a rheogram 
characteristic of pseudoplastic fluids. For n > 1, a curve characteristic of dilatant fluids is 
found. The index n, thus, is a measure of the degree of non-Newtonian behaviour of the 
fluid and may often be regarded as constant over several decades of shear rate. Now, for a 
Newtonian fluid, the viscosity is given by Equation (4-2). Analogously, one can define an 


"apparent viscosity" for power-law fluids by combining Equations (4-2) and (4-4) to obtain 
Boga Nahas (4-5) 


The power-law model is an attempt at empirical curve-fitting with maximum 
simplicity. Even though Equation (4-4) may fail to fit the total range of possible shear rates 
for some fluids, the expression can be very useful for a two-parameter fit of rheological data 


over a wide range of shear rates. One of the main objections to the power-law model is that 
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it does not predict the limiting apparent viscosities (i.e., the apparent viscosities at zero and 


infinite shear rates). 
4.1.3 Generalization of Darcy's Law 


The theory of laminar flow of Newtonian fluids through homogeneous porous media 
is based on the classical experiment of Darcy. Using the modified Blake-Kozeny equation 
for one-dimensional flow of power-law fluids through porous media (Christopher and 
Middleman, 1965; Gaitonde and Middleman, 1967; Savins, 1969), the superficial flow 
velocity can be expressed as (Ikoku and Ramey, 1979; Pascal and Pascal, 1985; Torok and 


Advani, 1987; Olarewaju, 1992) 


1/n 
ig eae (4-6) 
Hef L 


where the "effective viscosity” (Ug) is given by (Christopher and Middleman, 1965; 


Gaitonde and Middleman, 1967; Savins, 1969; Ikoku and Ramey, 1979; Olarewaju, 1992) 
Jaf n (1-n)/2 47 
Mest = 75 (9+3/nJ" (150k9) (4-7) 


It is clear from Equation (4-7) that in the Newtonian limit (1.¢., 2 = 1) Meg = H = Newtonian 
(constant) viscosity. From Equation (4-6), an analogy to Darcy's law for power-law fluids, 
neglecting gravity, can be expressed as (Ikoku and Ramey, 1979; Pascal and Pascal, 1985; 
Olarewaju, 1992) 


eapahoee, CP. 


= (4-8) 
mi Leg or 


where u, is the superficial velocity in the horizontal (radial) direction. 
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4.2 Flow Properties of Fractal Fracture Networks: Theoretical Background 


When a naturally fractured reservoir is highly disordered and fractal, the geometric and 
transport properties of the fracture networks differ in a non-trivial way from those for the 
corresponding Euclidean flow media (Chang and Yortsos, 1990; Acufia er al., 1992a; 
Acufia et al., 1992b). The theoretical, ideal response of perfect fractal objects is described as 
follows: many of their basic properties, defined as averages over a region of scale r, are 
scale-dependent and are proportional to non-zero powers of r. For example, the mass 
density of an arbitrary fractal network of fractures around an arbitrary point decreases in a 
power-law fashion with respect to the distance r (Acufia et al., 1992a; Acufia et al., 1992b). 
The exponent of this power-law is (d;— d), where dris the "mass fractal dimension" of the 
fractal network of fractures and d is the Euclidean dimension of the medium in which the 
fractal object is embedded and is an integer (1, 2 or 3). The physical meaning of the "mass 


fractal dimension", in the context of well testing, will be discussed subsequently. 


The calculation of fracture density is considered in the various fracture networks 
shown in Figure 4-1 (Acufia et al., 1992b), with the embedding medium being two- 
dimensional (so that d = 2) in every case. The concentric circles are at various radii from the 
centre (where the well is located) and the fracture densities are calculated within these 
arbitrary circles. The first network (network (a)) is nothing but a single fracture, the mass of 
which increases in direct proportion to the distance r, with the area (volume) of the 
embedding medium increasing in proportion to r¢. Thus, the density of the fracture (or the 
permeable flow path) within the two-dimensional medium varies in proportion to r ae Olean 
more general terms, to r4f—d where df= 1 and d = 2. Figure 4-1(c) exhibits a fracture- 
matrix system similar to the Warren and Root type double-porosity model with the 
horizontal and vertical lines representing the uniform fracture system. In this case, the 


fracture network is equivalent to a two-dimensional homogeneous medium of Euclidean 
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(c) dp =2 


Figure 4 — 1: Three Different Fracture 
Networks Embedded into a Medium of 
Dimension d = 2 (After Acuna ef al., 1992 b) 
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geometry (Chang and Yortsos, 1990). For this network, the fracture density, again 
proportional to r4f~ 4, is invariant with distance because of the fact that dg=d=2.In 
network (b), the indicated mass of dotted lines is assumed to represent a fractal fracture 
network. This network corresponds to a non-Euclidean case where, even though the same 
power-law variation of fracture density is exhibited, the value of dr lies between one and 
two; in other words, in this particular case, the cumulative "mass" or length of fractures 
contained within an observed area grows faster than it does in the case of a line (a single 
fracture; network (a)) but not as fast as in the case of a plane (a Warren-Root type 
homogeneous fracture system; network (c)). In a three-dimensional flow space, a value of 
flow dimensionality (dp) less than or greater than two depends on a combination of wellbore 
geometry and the nature of the fracture network (Acufia et al., 1992a; Acufia et al., 1992b). 
For example, the fractal dimension may be less than two if the fracture network exists only 
in the horizontal plane and the formation is homogeneous in the transverse direction (Acufia 
et al., 1992b). Similarly, dgmay be greater than two if flow also occurs parallel to the 
wellbore as in the case of a partially penetrating well in a relatively thick reservoir (Acufia et 
al., 1992b). The parameter dris strictly a geometric property of a fractal object and can be 
estimated by measuring the values of cumulative mass (or length) of fractures contained 
within circular regions of known radii r (Acufia et al., 1992a). In fact, as is evident from the 
discussion presented in this paragraph, the slope of a linear log-log plot of cumulative mass 


(or length) of fractures against radius would be the mass fractal dimension, dy. 


Another parameter (of fractal objects) of interest is the fractal exponent, 0, which is 
related to the topology of the fracture network (@ 2 0) (Beier, 1990a; Beier, 1990b). This 
parameter characterizes both the geometric and transport properties of the fracture network 
(Acufia et al., 1992a) and, hence, has a significant role to play in the diffusion of fluid along 


paths constrained to fractal geometry. During diffusion in fractal objects, the mean square 
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distance (from the origin) of a single particle diffusing from the origin at time t is given by 


the following scaling relation (Orbach, 1986; Acufia et al., 1992a; Acufia et al., 1992b) 


ENE TENE, (4-9) 


When 6 = 0, Equation (4-9) reduces to <r*(t)> « t, which is the corresponding scaling 
relation for the "normal situation" of diffusion in Euclidean space (e.g., the cases shown by 
Figures 4-1a and 4-1c). In these two cases, the fracture connectivity is high and therefore, 
diffusion is easier. On the other hand, in the fractal network of fractures (represented by 
Figure 4-1b), there is a diffusion slowdown (for @ > 0) because of poor spatial fracture 
connectivity and because of a higher degree of tortousity of the fluid flow path (Orbach, 
1986; Beier, 1990a; Beier, 1990b; Acufia er al., 1992a; Acufia et al., 1992b). In general, 0 
increases with poorer fracture connectivity and higher tortuosity of the flow paths (Acufia et 


al., 1992b). 


The spectral dimension of a fractal object, d, (Beier, 1990a; Sahimi and Yortsos, 1990; 
Beier, 1990b), is an extremely important parameter of the flow network, especially in the 
context of pressure transient analysis. The parameters d; and @ are related to the spectral 


dimension d, in the following way (Orbach, 1986; Beier, 1990a; Beier, 1990b) 


d, = f (4-10) 


Beier (1990a) and Acuiia et al. (1992a, 1992b) have shown that the asymptotic slope of the 
transient pressure curve in a fractal reservoir is related to the spectral dimension d, and 
hence, in Beier's studies (1990a, 1990b), the parameter d, is used to describe the fractal 
network instead of @. For a totally areal flow network, d, is less than or equal to two (Beier, 
1990a; Beier, 1990b; Acufia et al., 1992b); such a situation may occur when the flow 


medium is homogeneous in the vertical direction or when the producing unit is relatively 
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thin (Acufia et al., 1992b). The single-well pressure transient behaviour during the flow of a 
single-phase fluid through a fractal reservoir may differ significantly from its homogeneous 
radial flow counterpart. The theoretical formulation of the pressure transient behaviour for a 
fractal fracture network was presented by Chang and Yortsos (1990) who used the 


following relations for fracture porosity ("mass density") and permeability 


d,—d 
of2) (4-11) 
0 


(4-12) 


o(r) 


and 


k(r) 
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where gg and kg are the porosity and permeability values at r = rp. The defining equation for 
the permeability of a fractal object was obtained in the study of Chang and Yortsos (1990) 
by using analogies from the work of physicists who studied the variation of conductance in 
fractal systems (O'Shaughnessy and Procaccia, 1985; Orbach, 1986). Similar definitions of 
porosity and permeability were also used in many other studies on pressure transient 
analysis for fractal reservoirs (Beier, 1990a; Beier, 1990b; Acufia et al., 1992a; Acufia et al., 
1992b). It is to be noted that these definitions of porosity and permeability are not point 
values, as traditionally understood, but macroscopic values over a region of size r. 
Obviously, in the Euclidean limit of dr = d and @ = 0, these porosity and permeability 
values become independent of scale and reduce to the corresponding point values. Equations 
(4-11) and (4-12) show that the conductivity and storativity of a fractal reservoir are power- 
law functions of scale with different power-law exponents. Thus, the diffusivity is also 
scale-dependent (in a power-law fashion) with a non-zero exponent which has been shown 
to be related to the topology of the fracture network (Chang and Yortsos, 1990; Beier, 


1990a; Sahimi and Yortsos, 1990; Beier, 1990b; Acufia et al., 1992a; Acufia et al., 1992b). 
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4.3 Assumptions 


The following assumptions are made in deriving the mathematical model considered 

in the present study: 

a) small pressure-gradients exist throughout the reservoir at all times, 

b) gravitational forces are negligible, 

c) the radial flow takes place through an areal fractal network of fractures with the 

well penetrating the entire formation thickness, 

d) the porous medium has a uniform thickness, 

e) fluid compressibility is small, and 


f) the non-Newtonian fluid obeys the Ostwald-de Waele power-law relationship. 


The scalings for porosity and permeability of the areally heterogeneous reservoir are 
given by Equations (4-11) and (4-12), respectively. From these equations, with rg = ry (the 
wellbore radius) and for d = 2, the spatial variations of porosity and permeability of a two- 


dimensional areal fractal network can also be expressed as: 


Ourir (4-13) 


$(r) 
and 


k(r) = kor ir, poe? (4-14) 


Equations (4-13) and (4-14) are similar to Beier's Equations (6) and (7) (Beier, 1990a), 
respectively. The only difference between these two sets of expressions is that in Equations 
(4-13) and (4-14), the reference length scale is the wellbore radius instead of the equivalent 


wellbore radius, as was used by Beier. 
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Now, from Equation (4-8), the superficial flow velocity for a power-law fluid can be 


expressed in radial coordinates as: 


Z kr) a I/n 
u, = {ae 7 2| (4-15) 


where [U.., (7) can be expressed, by combining Equations (4-7), (4-11) and (4-12), as 


(n-1)(d; /d,—d, +1) 
(4-16) 


Hig(t) = A (= 
: 


w 


where 
H eal 
A= 99 (9 + 3! n)" {150k bf 2 (4-17) 


It is to be noted that Equation (4-15) is a modified version of the generalized form of 
Darcy's law as has been considered in some previous studies (e.g., Ikoku and Ramey, 1979; 
Pascal and Pascal, 1985; Olarewaju, 1992). However, in all these studies the power-law 
flow is considered to take place through a flow medium (or media) having no spatial 
variability of its transmissive or storage properties. In the present study, the generalized 
form of Darcy's law for a fractal object within which a power-law flow is taking place is 
extended. Additional assumptions made in obtaining a linearized partial differential equation 


are discussed later in this chapter. 
4.4 Some comments about the use of Equation (4-15) 


Yortsos (1991) extended results presented in the physics literature to a study of non- 


Newtonian power-law flow in percolation systems. The theory of percolation, originally 
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proposed more than thirty years ago and thus named because of an analogy with the process 
of a fluid percolating through a solid, continues to be one of the important tools for 
physicists in the study of the geometric and transport properties of porous media (Sahimi 
and Yortsos, 1990; Yortsos, 1991; Balberg er al., 1991; Berkowitz and Balberg, 1993). The 
theory has been developed extensively in the field of statistical physics and has been applied 
in a variety of problems including the conceptualization of geometrical properties and 
transport phenomena observed in disordered flow media (Sahimi and Yortsos, 1990; 
Yortsos, 1991; Balberg et al., 1991; Berkowitz and Balberg, 1993). Percolation theory 
provides "universal" scaling laws which determine the physical and geometrical 
characteristics of a system (Balberg er al., 1991; Berkowitz and Balberg, 1993). These laws 


are generally expressed in power-law forms of the type (Balberg et al., 1991) 
A « (v-v,)” (4-18) 


where A is a geometrical or physically observable quantity, v is the fractional volume of the 
conducting phase and v, is the critical (threshold) value for the onset of percolation (i.e., of 
system connectivity) and x is the exponent for quantity A and can be determined from theory 
and/or computer simulation and/or experiment. Scaling laws of the type shown above hold 
for v relatively close to v¢ (typically, for v < 2v, (Berkowitz and Balberg, 1993)). For the 
universal scaling laws such as Equation (4-18), the exponent x is dependent only on the 


system dimensionality, d. 


Making use of a scaling law similar to Equation (4-18), Sahimi and Yortsos (1990) 
demonstrated the following scaling relationship, with length, of effective hydraulic 


conductivity of a percolation system for Newtonian flows: 


Kecgi a (4-19) 


23 


afl. < ng? calcd ~ 


i: +t - tee ® 4 
joe) sound, Grew lien aaanlietae du derqere’s Sueur 


saya) Sah. 4d CC} aeelinfh Pees ae “ST: , ts Veg Be dart) aster? aw 


Cues) Sarr a 


Sus ont site tks Pein up i | 


ap Ds 
0 


; a ay Sy 
roo sepmg Taaearanacinge be valluSiaeeehei. ot wb 


7 1 


erty tenants, CEU CE on sbarcat Ott in m yredioll 
yamosg. torus oni ie sel igidr wheat gedlua? “leat 


fhe Soe kere tt teat alee <* GoenmeyeS 


» on ro ST 


wip te storey fofeothd Aes 20 eat eal? culyruisdy Vilas veg 10 ‘S3°TRS 
i) oust Bleep’ “i nee ay pique pinnate) lesetica eer ob 7 Bas 
eed ei and ripidls Spine Cra Ran Bee rege? ot nd tia (gh i 

patin Ointh ee ay a Toews! shtiao® Sipser aes TONE 8 CH ti 
10h k FF KOON orehee beak cohiOR tS cht ee fr or ' 
Av te " fu Fotaayay oF HiMaee wl - voahis> 32 Pat owal 


a = wt ; - 
ee eee | Meee ie Int alg asain am gal 
Sugita soem hes tara vant ioe yen 


where K = 2 (Sahimi and Yortsos, 1990) and v = 0.88 (Sahimi and Yortsos, 1990) (v = 
0.854 in the study of Berkowitz and Balberg, 1993) for a 3-D system. It should be noted 
that in the study of Sahimi and Yortsos (1990), 4: was used instead of Kas the conductivity 
scaling exponent. For a non-Newtonian power-law flow in a percolation system, the 


conductivity scaling has been shown to possess the form (Yortsos, 1991) 
Reset ly (4-20) 


where t(n), although not precisely known, may be defined approximately for 3-D 


percolation clusters as (Yortsos, 1991) 
t(n) = 1.76 + 0.24/n (4-21) 


To the best of our knowledge, an equivalent expression for t(n) for 2-D systems is not 
currently available. For an areal flow system, as is considered in the present study, the use of 
a proper flow equation would require the availability of a suitable approximating form for 


t(n) so that the conductivity scaling given by Equation (4-20) may be used. 


In the present study, use has been made of Equation (4-15) to describe the flow 
velocity of a power-law fluid in an arbitrary fractal medium. From Equations (4-14) through 


(4-16), the scaling for non-Newtonian power-law flow in a fractal has been found to be 


ta cents 4-22 
or ( ) 


Tw 


dp-l-de/d,)+(/n)(l—d Id, 
( . } f pld,)+Q/n)l—d¢ eae 


It should be noted that at this point, to our knowledge, rigorous theoretical and/or 
experimental results for power-law flow in fractal objects have not been presented in the 
literature, and hence, the validity of the scaling presented by Equation (4-22) is questionable. 


However, the present study deals only with areal flow networks, and because of the lack of 
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available scaling exponents for power-law flow in 2-D percolation systems, use could not be 
made of a hydraulic conductivity scaling law similar to that considered in the study of 
Yortsos (1991). The rationale for proceeding in our present analysis with Equation (4-22) 
also lies in the fact that this equation arises from a straightforward coupling, as it were, of 
the flow equations used in studies dealing with non-Newtonian power-law flow in 
homogeneous systems (i.e., Ikoku and Ramey (1979)) and those dealing with Newtonian 
flow in fractal systems (e.g., Beier (1990a)). And, as will be shown later, the resulting 
partial differential equation and its solutions consequently are direct generalizations of those 
presented in the aforementioned studies, and they also reduce to the classical pressure- 


transient models in the limit of n = 1 and dy = d,=2. 
4.5 Partial Differential Equation 


The continuity equation for radial flow in a porous medium may be written as 


Wh Ya) 0 
= 5, ("PH ) = 5 (PP) oe 


For isothermal fluid flow and assuming constant fluid compressibility, the density of the 


fluid at any pressure, p, can be expressed as 
P = Lp ec(P-Po) (4-24) 


where p, is the density at a given pressure, py. From Equation (4-24), one obtains 


eae C2 (4-25) 
or a or 

and 
Caaae.. oP (4-26) 
eames 3 
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Also, we have 


re) 
a af ob (4-27) 


where Cr is the effective pore space compressibility. Introducing Equations (4-15), (4-25), 


(4-26) and (4-27) into Equation (4-23), one gets, after simplification 


I/n I/n l+n 
LO) | _k(r) op k(r) 2B)" : 2 
£3 | io 4 |: | & ean Naa 


n—-] 


Multiplying both sides of Equation (4-28) by [xr 2) " | it is possible to obtain, after 
ie 


further simplification 
: op op ad -lin 1 
malt © | EAT) or ar’ ose np! (r a (r) ze 
n-1 
k(r) = 2 dp 
Meg ( arin 4 P ered ke ar} ot (4-29) 


If constant and small fluid compressibility and small pressure gradients are assumed, the 


gradient-squared term may be considered negligible and Equation (4-29) reduces to 


18) B] + [ EAE 
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n or i idee ee t| 27. 20 
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a Op| » op 
= O(r)cfLin Ani 2) Or (4-30) 
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Equation (4-30) is a nonlinear partial differential equation and represents the governing 
equation for the radial flow of a non-Newtonian power-law fluid through an areally 


heterogeneous porous medium with a fractal structure. 


For a Newtonian fluid, n=1 and w,. = H = pw, and for this case, Equation (4-30) 


reduces to 


dediler re OP ae 
Z 2 ret P| - = (rou 2 (4-31) 


which is the diffusion equation for radial flow of a slightly compressible Newtonian fluid 
through a heterogeneous reservoir (e.g., Equation (5) of Beier, 1990a). Again, by setting dy 


= d, =2(so that, k(r) =k, o(r) = 0 and Leer) = 1), Equation (4-30) becomes 


n—-1 
0’ p nop Log ) (2) n Op 
ad din ied Sli ae al Shoal cd 4-32 
a3 nde k or Ot pee) 


which is the diffusion equation for radial flow of a non-Newtonian power-law fluid through 


a homogeneous porous medium (similar to Equation (A-9) of Ikoku and Ramey, 1979). 


Equation (4-30) is linearized by making an approximation suggested by Ikoku (1978) 
and Ikoku and Ramey (1979). From Equation (4-15), one has 
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where q is the production (or injection) rate. Substituting Equation (4-34) into Equation (4- 
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The approximation, given by Equation (4-34), has allowed us to linearize Equation (4-30) 
so that proper (albeit approximate) analytical solutions can be obtained. This approximation 
is equivalent to assuming that the flow rate at any radial distance from the wellbore is a 
constant. For non-Newtonian power-law flow in a homogeneous porous medium, Ikoku 
(1978) and Ikoku and Ramey (1979, 1982) have shown that the analytical solution obtained 
by making such an approximation compared fairly well with a more rigorous numerical 
solution when n is between 0.5 and 1. In the present study, it has been observed that for 
power-law flow through fractal reservoirs, such an approximation results in errors that are 
not large for many values of d, and n. A detailed discussion of this aspect of the analytical 


model will be presented in a later chapter. 


Now, introducing Equations (4-13), (4-14) and (4-16) into (4-35) and simplifying, 


one obtains 


where 


a mAs (2) (4-37) 
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4.6 Dimensionless Variables 


The dimensionless groups are defined in the following manner: 


a (p;— Pp) Kk ee d,/d, 


a A (q/2mh)" Grek 

es 4-39 

ee (4-39) 
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It should be noted that the dimensionless groups defined in Equations (4-38) through (4-40) 
are quite similar to the corresponding groups defined by Beier (1990a). In fact, for n = 1 (so 
that A = H = p, and G = pg,,c,/k,,), the two sets of groups are identical. Moreover, for dy = 
d, = 2 (so that A = Leg) Equations (4-38) through (4-40) reduce to the corresponding 
dimensionless variables, defined by Ikoku (1978) and Ikoku and Ramey (1979, 1980, 1982) 


for the flow of a power-law fluid through a homogeneous porous medium. 


From Equations (4-38) through (4-40), it can be shown that 
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Substituting Equations (4-41) through (4-43) into Equation (4-36) and simplifying, one gets 
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(4-44) 


Equation (4-44) is the dimensionless form of the linearized partial differential equation for 
the flow of a power-law fluid through a fractal reservoir. For a homogeneous porous 
medium (dy = d, = 2), Equation (4-44) reduces to the dimensionless linearized partial 
differential equation (Equation (B-1) of Ikoku and Ramey, 1979) for the flow of a power- 
law fluid through a reservoir. Also, for n =1 (Newtonian fluid), Equation (4-44) is the same 


as Beier's Equation (A-1) (1990a). 
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Chapter V 


SOLUTIONS FOR THE INFINITE RESERVOIR CASE 


The emphasis in the present chapter will be on discussion of the solution of the linear 
partial differential equation (Equation (4-44)) for the case of constant production rate from 
an infinite reservoir. Wellbore storage and skin effects will not be considered in the present 
analysis. Subsequently, the rate solution of Equation (4-44) for the case of constant wellbore 


pressure in an infinite system will be discussed briefly. 
5.1 Constant Rate Case 
5.1.1 Initial Value Problem: Finite Radius Well 


The case of production of a power-law fluid at a constant rate from an infinite reservoir 
into a finite-radius wellbore is considered in this section. The initial value problem becomes 


(Equation (4-36)) 


2 d A {(dy/d, (nt1)-d, (n-1)+n-3} 
a a ip eae n op = ca <i op (5-1) 
or Tk aOR hel ¢ absicad as ot 
with the initial condition given by 
p(r,0) = p,,forallr (5-2) 


The inner boundary condition is given by the following (from Equations (4-6) and (4-7)) 


relationship 


op 
or 


Pee (5-3) 
r=Try ky, 2mr,,h 
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The reservoir is infinitely large so that the pressure at any time and at all locations far away 
from the production well remains essentially constant at the initial pressure, p;. Thus, the 


outer boundary condition is 


lim p(r,t) = p,, forallt (5-4) 


In terms of dimensionless variables, Equations (5-1) through (5-4) become, respectively 


Z 
Oo Pp Ji. E Sar aera = (a, 14, ) rp teste nt-aytnnt-2} Pp 


Or,” meer nd. dr or. Oty 
(O=)) 

Dito) = 0, forall rp (5-6) 

OP p d, 

— = ——,forallt eo 

ro |, a. D (5-7) 

and 
lim Py(% tp) = 0, forall tp (5-8) 


Equation (5-5) is a linear parabolic partial differential equation expressed in terms of 
dimensionless variables and may be solved by applying the Laplace transformation. 
Application of the transformation to Equation (5-5) and the initial and boundary conditions 


yields 


25 d d n+i})— y= 
pe Pp ty c 4), Be ce (d, 1,) ry" 1)=d,( hp 


(5-9) 
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and 


lim Bp(tp.1) = 0 (5-11) 


with the initial condition having been used in obtaining the transformation of the time- 
derivative of pp in Equation (5-9). By a suitable change of variables, Equation (5-9) can be 


transformed into a Bessel equation. Let us assume that 

Pp = f(p) (5-12) 
where 

Sas (5-13) 


and where a@, B and y are arbitrary quantities to be determined. From Equations (5-12) and 


(5-13), it can be shown that 


hy =. ig Sagage (5-14) 
dr, is i 

and that 
c= 2 ny 
ds EE ea Bg rs2yyiop + Me ae (5-15) 
ary «< es an 


Substituting Equations (5-12) through (5-15) into Equation (5-9), one obtains, after 


simplification 
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(5-16) 


It is now possible to choose values of the parameters a, B and y such that Equation (5-16) 


reduces to a form for which standard solutions are available. Choosing 


2V1(d, /d_) 
= og ee ene (5217) 
(d,/d,)(n+1)-—d,(n—-1) 
ad-id +1])-d,(n-I1 
g = (Ala jnt1-d)(n=) ae 
Z 
and 
d./a.)(n+1)—nd 
Y= (A nae (5-19) 
2 
it can be shown that Equation (5-16) reduces to 
2 
ee | fei (5-20) 
n+I1-d,(n—I) 
Equation (5-20) is Bessel's modified differential equation of order 
v=]- Sa (5-21) 
n+1—d,(n-1) 


The magnitude of the parameter v is quite important in the present analysis. As will be 


shown later, the value of v represents the large-time slope of the log-log plot of 
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dimensionless wellbore pressure against dimensionless time for the flow system under 


consideration. 


The general solution to Equation (5-17) is (Abramowitz and Stegun, 1972) 
f(p) = Cl,(p) + CK,(p) (5-22) 


where /,, and K,, are modified Bessel functions of the first and second kind, respectively, and 


of order v; the parameter v can have integral or nonintegral values. By combining Equations 
(5-11) and (5-22), it can be shown that C; = 0. Thus, from Equations (5-12), (5-13) and (5- 


22), one can write 


Baca C7) K (ars) (5-23) 
so that 

Po = C, [ yK, (a) + a@BK’, (a)] (5-24) 

ary Be 


Using the following property of the modified Bessel function, K,, (Carslaw and Jaeger, 


1959): 
K',(@) = -—K,(0)-K,.,(@) (5-25) 


and noting from Equations (5-18), (5-19) and (5-21) that Bv = y, Equation (5-24) can be 


rewritten as 


4) K)_,(@) (5-26) 
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The constant C2 can be determined by comparing Equation (5-26) with the inner boundary 


condition (Equation (5-10)) and is given by 


1 


C= => 
; IZ Keio) 


(5-28) 


Thus, from Equations (5-23) and (5-28), one can write the solution for Equation (5-9) as 


PaKE( OTs ) 


pe K 


Pare t ee 
p(T,!) Ae) 


(5-29) 


where v is defined by Equation (5-21), and a, B and y are given by Equations (5-17), (5- 
18) and (5-19), respectively. Equation (5-29) is the Laplace transform of the general solution 
for the transient pressure distribution during the constant rate flow of a power-law fluid in an 


infinite reservoir showing fractal behaviour. 


The general solution at the wellbore (rp = 1) is given by 


Ky(@) 
pee K,_,(@) 


cj 

n+I1-d,(n-1) 

Eas 2v1 
*|n+1-d,(n-—1) 


ey) 


(5-30) 
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Beier (1990a) concluded that the dimensionless wellbore pressure during the production of a 
Newtonian fluid from a fractal reservoir depends on d, and is independent of dy because of 
the way the dimensionless variables were defined in his study. From Equation (5-30) it can 
be seen that the dimensionless wellbore pressure during the production of a power-law fluid 
from a fractal reservoir depends on d, and n and is independent of dy. Part of the dependence 
on drhas been absorbed in the definition of the dimensionless groups. The linearization 


approximation is also responsible for the dimensionless wellbore pressure solution being 


independent of dp. 


For zero wellbore storage and Newtonian flow (n = 1), Equation (5-30) reduces to 
Equation (A-13) derived in Beier's study (1990a). Also, for d= d, = 2, Equation (5-30) 
exhibits the wellbore pressure response during the flow of a power-law fluid through an 
infinite homogeneous porous medium (Equation (B-13) of Ikoku and Ramey, 1979). 
Equation (5-30) may be evaluated directly by numerical inversion. However, it is also 


possible to find approximate analytical inversions, as will be shown later in this chapter. 
5.1.2 Limiting Solutions for Early and Late Times 


It is possible to derive the early- and late-time approximating forms of Equation (5- 
30). The limiting form for small times will be presented first. For early times, one has 


1 —> oo , when the modified Bessel function, K,(z), can be approximated as (Abramowitz 


and Stegun, 1972) 
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where z | = ————————__ 
n+I1-d,(n—1) 


— co . Substituting Equation (5-31) into Equation (5-30) and 
simplifying, it follows that for ] — oo , Equation (5-30) can be approximated as 


Dull) eae (5-32) 


Inverting Equation (5-32) analytically (Abramowitz and Stegun, 1972), it can be shown that, 


at early times, the dimensionless wellbore pressure can be expressed as 


~ 


Pyp(tp) = 2 |-> (5-33) 
1 
Equation (5-33) also represents the early-time pressure solution for a finite-radius well in a 


homogeneous formation (Streltsova, 1988). This early-time pressure response, when plotted 


against time on logarithmic coordinates, exhibits a straight line with a slope of 1/2. 


Similarly, one can also consider the dimensionless wellbore pressure response (given 
by Equation (5-30)) at large times. An expression for the modified Bessel function, K,(z), 


can be written as (Abramowitz and Stegun, 1972) 


a lez) "i. DG) 


a ae 2 sin (vv) 
Z ne [1.,(z) — 1,(2)| (5-34) 


where I'(a) is Gamma function and where (Ikoku, 1978) 
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Thus, from Equations (5-34) through (5-36), one gets 


K,(z) = eyez) ie: (2/2) | ‘, 
i ? 2 ee wer ~ 


el i a (zfey v2 
2 Vil+v) F(2+v) 0 (5-37) 


For large times, such that / (and hence, z) is very small, Equation (5-37) can be 


approximated by neglecting z2 and higher powers of z within the winged brackets to obtain 


K.(z) = soe ey i ie) i 
. 2 2 iio) 2 Ly) 


_ roe)” re (z)" T'(1-v) (5-38) 
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From Equation (5-38), it is possible write 
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Neglecting z? and higher powers of z within the square brackets, Equation (5-39) reduces, 


for small magnitudes of v, to 


Ky (2z) =e I'(v) (2) 2 q T(v) a 
Koz) el (l=v)\2 a aes ESC a (5-40) 


Applying Equation (5-40) to Equation (5-30), one gets the large-time approximation of the 


wellbore pressure solution, in Laplace space, as 


eo = d, (n+1)(d,—1) 
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(i) 
ered = 
I(n+I—nd,) (5-41) 


Inverting Equation (5-41) term by term using Laplace transform tables (Abramowitz and 


Stegun, 1972), one gets 
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Bed = C(t 1) 


At sufficiently large times, one can obtain an approximate expression for p,,n, from 


Equation (5-42), as follows: 
(nt1)(d,-1) 
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n+1l—nd 
Pe ] 
(nt+1—nd,) (n+1-nd,) (5-43) 


Equation (5-43) may be compared with dimensionless wellbore pressure values obtained by 


numerical inversion of Equation (5-30). 


Equation (5-43) shows that for practical purposes (i.e., at relatively large times) a 
n+1—nd, 


n+1—(n—1)d, 


Cartesian graph of p,,p against tp will yield a straight line with a slope, m, given by: 


(n+1)(d,-1) 


iy n+1—d,(n-1) | 1 nt1-d,(n—1) ] 


Tt —————— (5-44) 
n+l—d,(n—I) n+d—nd, 
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and an intercept, at tp = 0, of -I/(n+1-nd,). 
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Another late-time (/ — 0) approximating form of Equation (5-30) may be derived. When 
z— 0, the modified Bessel function, K,(z), can be approximated as (Abramowitz and 


Stegun, 1972) 


il pane 
K,(z) = Sr(vi(2) (5-45) 


Using Equation (5-45), Equation (5-30) can be reduced to the following form: 
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a i ia) 
(5-46) 
Inverting Equation (5-46) by using Laplace transform tables (Abramowitz and Stegun, 


1972), one obtains, after simplifying, 


(n+1)(d,-1) 
Renee (Ted Hl nes dined 
Pwo (t5) ~ _ntl~d(n-1) | SY peters Se 
- d. eee 
n+1—d,(n—1) 
(n+1-nd, ) 
pee (n—-1) 
15 11 a ata (5-47) 
nei —nd, 
: : : tl : 
Equation (5-47) is the same as Equation (5-43), except that the last term, Sam remy" 1s 


omitted from the former equation. Equation (5-47) would apply, therefore, at very large 
times and for positive values of the exponent of dimensionless time. It can be noted that the 
exponent of tp in Equations (5-43) and (5-47) is given by the parameter v (see Equation (5- 
21)). Thus, for positive values of v and at large times, a log-log straight line plot of pyp 


versus tp would exhibit a slope of v. For a pseudoplastic fluid (0 < n < 1) and DS eeaAy 
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would always be positive. However, for a dilatant fluid (1 < n < 2), in order for v to be 


poSitive, it is necessary that d, < 1.5. 


It is obvious that for n = 1 (Newtonian fluid), a large-time plot of the pressure 
transient, given by Equation (5-47), on logarithmic coordinates would yield a straight line 
with a slope of 1 — d;/2, from which the value of d, can be obtained. A similar conclusion 
was also made by Chang and Yortsos (1990) and Beier (1990a). It has also been shown that 
at large times an identical behaviour would be exhibited by a line-source well in a fractal 
reservoir (Acufia et al., 1990a, 1990b). Equation (5-47) also shows that for a homogeneous 
reservoir (df = ds = 2) and for a pseudoplastic fluid, a large-time log-log plot of pp versus 
tp would be characterized by a straight line with a slope of (1—-n)/(3-n), as has been 


observed earlier by Ikoku and Ramey (1978). 
5.1.3 Pressure Solutions for Line-Source Wellbore Case 


It is possible to obtain transient solutions for dimensionless pressure not only in 
Laplace space but also in real space for the case of a line-source (sink) wellbore in an infinite 
fractal reservoir. In every respect, other than the inner boundary condition, this case is 
mathematically the same as the finite-wellbore case. For the infinitesimal-source case, the 
inner boundary condition may be described in dimensionless terms, from Equation (4-15), 


as 


yr +1-(dy/d, )(n+1) OD» 
D 


= -d,/d, (5-48) 


Ory rp 0 


Now, for the infinite outer boundary condition, the solution to Equation (5-9) is given by 


Equation (5-23) which includes the constant C2 to be determined from the inner boundary 


condition. From Equation (5-23), it can be shown that 
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pds td 1d (ntl) aD» 
D a 


dr, 


= Cr’ apK,(orp)| (5-49) 
Tp 0 3 


Using the following relationship (Barker, 1988) 


lim 2°K,(z) = 2 iy) fora 0 (5-50) 


and combining it with Equation (5-48) expressed in Laplace space, one can obtain an 
expression for Cy from Equation (5-49). The resulting expression for dimensionless 
pressure is given by 


Pp(%,!) — 


d Y ! B 
£2 (2) Al) yc (5-51) 


d, pBY(1-v) la Lee 
where a, 8, yand v are given by Equations (5-17), (5-18), (5-19) and (5-21), respectively. 
Equation (5-51) defines the pressure response (in Laplace space) of an infinite fractal 
reservoir due to a constant rate of production from it by a line-source well. Equation (5-51) 
can be inverted to the real plane analytically; this analytical expression can be written as 


follows 


27 Doe} 
joni Oeplin) aa 5 r-» aa (5-52) 


where the incomplete Gamma function is defined as 
Piyap ety ex ea (5-53) 


and the parameter a = d7/2 Bd,. At the wellbore, Equation (5-52) reduces to 
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Pao(tp) = reer r-» =| (5-54) 


tp 


It can be seen from Equation (5-54) that the wellbore pressure response at a given time 
depends only on n and d,, as in the finite-wellbore case. The incomplete Gamma function 


can be expressed as follows (Barker, 1988) 


-I 
{b,x} = T(b) - ee en S (5-55) 


For small values of the argument a2/tp, which occur after a short time, one may use only the 
first two terms of the expanded form of the incomplete Gamma function (Acufia et al., 


1992b) and thus, from Equations (5-54) and (5-55), it can be shown that 


ai” 7 


oO + 
ei) 7 bpd (ne) 


(5-56) 


Equation (5-56) is the large-time approximation for the line-source wellbore pressure 
response. A similar equation has been derived also for the case of n = 1 by Acufia et al. 
(1992b). It is interesting to note that Equation (5-56) is, in fact, identical to the large-time 


approximation for the finite-wellbore pressure solution (given by Equation (5-43)). 


For n = 1 and dy= d, = 2, Equation (5-52) reduces to the following form 


aalectid) Ste E(B (5-57) 
Zaret,, 


where the exponential integral (Theis well function) is given by 


E(x) = i du (5-58) 
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5.1.4 Skin Factor 


From Equations (4-38), (4-40) and (5-43), one can write 


DRED aD): 
2v : 
: 1 (4, /d,)(n+1)-d,(n—D)} t fee 
C,(n+1—nd,) (Gres ity) , ee 
where 
rok (d id 
~ Se (5-60) 


PP A(q/2ah)" 
and G is given by Equation (4-37). 


The van Everdingen and Hurst skin factor is defined as a dimensionless constant, s, 
which relates the pressure drop in the skin to the dimensionless rate of flow. The skin 


pressure drop can be defined (Ikoku, 1978) as 


Ap, = — (5-61) 
C, 
Combining Equations (5-59) and (5-61), the total pressure drop can be expressed as 


Ap = D; — Pug 


2v 
1 | {(4,/4,in+-a(n-1} 0” elae 
C,| (nti-nd, Gri") T(1-v) — (n+1=nd,) 


+5 (5-62) 


Att=0, Ap = Ap, and the skin factor can be calculated from the above equation as 
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s = Ap,.C, + ———_ -63 
AOS: n+1-nd, m2) 


It is to be noted that, for a fractal reservoir, the value of the skin factor obtained (by using 
Equation (5-63)) would also include any near-wellbore manifestations of deviations from 


the fractal distributions of the reservoir properties. 
5.1.5 Radius of Investigation 


Under steady-state conditions, Equation (5-5) reduces to 


— = 5-64 
dr,” ok Males dole Coney 
Integrating Equation (5-64) once with respect to rp, one obtains 
meres | Ad 
: J 2 = 5-65 
is A dr 1 ( ) 
where Aj is a constant to be determined. The boundary conditions are 
A Sh ate (5-66) 
dr, eT d, 
and 
Pp = 9@% = Vi (5-67) 
where rp; = ‘iny/y. From Equations (5-65) and (5-66), it can be shown that 
ies, lao (5-68) 
id d, 
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Introducing Equation (5-68) into Equation (5-65) and integrating once more with respect to 


rp, one obtains 


ACE Id, )(n+1)—nd, 
pf cetera aoe 


Pra Ay (5-69) 


n+I1-nd, 


where A> is a constant to be determined. Using the second boundary condition (Equation (5- 


67)) and Equation (5-69), one can determine A> to be 


(dy ld, )(n+1)—nd, 


i ees ons 5-70 
; n+1-nd, a 


Combining Equations (5-69) and (5-70), the dimensionless wellbore pressure can be 


expressed as 


iy errs ¥ 1] 


ee 5-71 
Pwo n+ ] z nd, Di ( ) 
Comparing Equation (5-71) with Equation (5-43), one obtains 
(d,/d,)(n+1)—nd {n Fied(n— py tp 
ias'4, ae ne REE oe Woe (5-72) 


T(1-v) 
from which the value of the radius of investigation at a given time can be calculated. 


For dr = d, =2, Equation (5-72) reduces to 


2 1/(n-1) 
on ieee ae irs = -)| (5-73) 


so that 
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where G is defined by Equation (4-37). Equation (5-74) is the expression for the radius of 
investigation derived by Ikoku and Ramey (1979) for power-law flow in a homogeneous 


porous medium. Moreover, for n = 1, Equation (5-74) becomes 


rauuup ny mlat (5-75) 


which is the radius of investigation of a constant-rate test during the flow of a Newtonian 
fluid of viscosity H through a homogeneous porous medium of permeability k,, and 


porosity 0, 
5.2 Constant Pressure Case 
5.2.1 Initial Value Problem 


The linear partial differential equation governing the transient flow of a non-Newtonian 
power-law fluid through a fractal reservoir is governed by Equation (5-1) expressed in 
terms of dimensional quantities. The initial and outer boundary (infinite) conditions are 
given by Equations (5-2) and (5-4), respectively. For a constant-pressure inner boundary, 


the inner boundary condition is expressed as 
Dir =7,, 1) = p, a) 


The dimensionless radius and time are defined by Equations (4-39) and (4-40), respectively. 


For the constant-pressure inner boundary condition, the dimensionless pressure is defined as 
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D; re Be 
It can be shown easily that using this new definition of dimensionless pressure and 
Equations (4-39) and (4-40), Equation (5-1) can be rewritten in dimensionless form as 
Equation (5-5). In order to solve Equation (5-5), one applies the Laplace transformation not 


only to the equation but also to the relevant boundary conditions. For the constant-pressure 


wellbore case, the inner boundary condition can be expressed in Laplace space as 


a 1 
Pp (% = 1,1) = a (5-78) 


As has been shown for the constant-rate case, the equation for dimensionless pressure (in 
Laplace space) after applying the initial and the outer boundary conditions is of the form 
given by Equation (5-23). Finally, by using the inner boundary condition, the pressure 
solution is obtained as 


re K, (ars) 


eee OPK: (a) 


(5-79) 


Let us define, for the constant-pressure inner boundary case, the dimensionless rate in 


a manner similar to that of Poon and Kisman (1992); that is 


Ar q : 
a, oll 5-80) 
2 Rep (D; —D,,) d; nO: =] : 


where the dimensionless rate is related to the pressure gradient at the wellbore as 


qo(d,/d,) = - ae, 


5-81 
Or, ( 


Tp=l 
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Taking the Laplace transform of Equation (5-81) and using Equation (5-79), the following 


expression for the dimensionless rate is obtained 


= BOGS. 5 
a MT K,(a) (5-82) 


Defining the cumulative production (over a given time tp) as 
Oy = JQ aty (5-83) 


and taking its Laplace transform, one gets, from Equations (5-82) and (5-83), the following 


expression for the cumulative production in Laplace space 


toy | ar 5-84 
2p Wl K,(a) Os 
5.2.2 Limiting Solutions for Early and Late Times 


At early times (/ — oe ), the approximation of the modified Bessel function, K,(z), is given 
by Equation (5-31). Substituting Equation (5-31) into Equation (5-82) and simplifying, it 


follows that for ] 4 ce , Equation (5-82) can be approximated in the real plane as 


| 


dp = xt, 


(5-85) 
Equation (5-85) represents the early-time approximation of the dimensionless production 


rate. Similarly, from Equations (5-31) and (5-84), the early-time approximate form of the 


cumulative production (in real space) can be obtained as 


Oh = 2,| (5-86) 
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The late-time (/ — 0) approximating forms of Equations (5-82) and (5-84) may also 
be derived. Using the approximation of the modified Bessel function, K,(z), when z > 0, as 
given by Equation (5-45), one can obtain the late-time limiting forms of dimensionless 


production rate and cumulative production, respectively, as 


’ oo 2v-1 faa? 
ape d. T(v) (5-87) 
and 
Tesvoetilens vt,” 
ae a ee De z 
Op = d, (I1—v)T(v) Pty, 
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Chapter VI 


NUMERICAL PRESSURE SOLUTION OF THE NONLINEAR EQUATION 
GOVERNING POWER-LAW FLOW THROUGH FRACTAL POROUS MEDIA 


In Chapter IV, a nonlinear partial differential equation (Equation (4-30)), governing the 
radial flow of a power-law fluid through a fractal network of fractures embedded into a two- 
dimensional medium, was derived. A "linearizing approximation" reduced this equation to a 
linear but approximate form (Equation (4-36)) and an analytical solution of the latter 
equation was subsequently obtained. The approximation is equivalent to assuming that the 
flow rate is time-independent at each radial location. Such an assumption is nearly correct in 
an expanding/contracting region close to the wellbore, but is not correct near the radius of 
investigation (Ikoku, 1978). This chapter presents a finite-difference scheme to solve the 
nonlinear partial differential equation for the flow of non-Newtonian power-law fluids 


through fractal reservoirs. 
6.1 Dimensionless Nonlinear Equation 


The dimensional form of the nonlinear equation is given by Equation (4-30). Making 
use of the definitions of the dimensionless variables (given by Equations (4-38) through (4- 


40)), and by using the following logarithmic transformation, 
x = In(rp), (6-1) 
Equation (4-30) may be written in dimensionless form as 


n-I 
OD, OD» | Bo.) * OPp : 
BED mee AoA = PD. | 4. SPD. (6-2) 
Sort ie tic sa an 
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where 


A, = n(d,-A,) (6-3) 
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(6-4) 
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(6-5) 


The initial condition and the inner (constant rate) and outer boundary (infinite reservoir) 
conditions are given, in terms of variables rp and tp, by Equations (5-6), (5-7) and (5-8), 


respectively. These equations can be rewritten, using Equation (6-1), respectively, as 


Dp(X, 0) = 0 (6-6) 
Pp = _ as (6-7) 
Ole, d, 

and 
lim pp(X, tp) = 0 (6-8) 


The numerical solution of the problem posed by Equations (6-2) and (6-6) through (6-8) is 


now inspected. 
6.2 Douglas-Jones Predictor-Corrector Method 


Ikoku (1978) and Ikoku and Ramey (1982) used the Douglas-Jones predictor- 
corrector method to obtain numerical solutions of the nonlinear partial differential equation 


governing the transient flow of a power-law fluid in a homogeneous reservoir. In this study, 
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use is made of the same method to solve Equations (6-2) through (6-8). Each of the finite- 
difference equations advances the solution by one-half of the time increment. In the 
predictor, the unknowns occur at the (j+ 1/2)th time level and the equations are linear. In the 
corrector, the unknowns occur at the (j+1)th time level and the equations are again linear. 
The predictor followed by the corrector is unconditionally stable; moreover, the systems of 
equations are of the tridiagonal form, and are easy to solve. The Douglas-Jones predictor- 


corrector method is second-order accurate, whereas the computing effort is doubled at each 


step. 


In order to solve a nonlinear partial differential equation of the form 


Oe ean aioe cB (6-9) 


= == ap I Uo [Dy 
ox? x Ot a( e ox 


use is made of the predictor (to advance the solution from the jth to the j+1/2th time level) 


given by 
EMDR a yg yf PnP | PE =P 
(Ax) i?" 7+1/29 Fi? 2Ax At/2 
j Bay 5 De (6-10) 
+ 8) Xs Ujeia» Pi» SAR 


and the corrector (to advance the solution from the j+1/2th to the j+1th time level) given by 


jee vi . Lites ; 
men —2p; aCe Diet ) t 5 (Pir — 2p; + Pis1) i j+i/2 Die me sin 
5 =f| Xstise Di 
(Ax) amt 
jel os j ' j+ll2 tee j+1/2 
en SS Er At Pi te i{» Di .1/2? vee Pist = (6-11) 
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Applying the predictor-corrector method to Equation (6-2), one obtains the predictor 


j+li2 j+ll2 jell2 j ; 1/2 
Di.) — 2D; Pe) = A,e™ te Ee = Pi iY (@ — p} | 


(Ax) 2Ax UNL 
+A, Fir Pi (6-12) 
and the corrector 
| jtl a Je fee et Jt+l 1 J 7 j 
3 =(p; -1 ~4P; + Diss ie 5 wa — 2p; + Diss) ah pean 
(Ax) 2Ax 


n-1 
eae pe Doe 5 pi’ - pi p! 
= ee A ONE? SEE 6-13 
Ao 2Ax At ee 


mete t — 1, 2, 3, ..., N+1 and j = 0, 1, 2, 3, ..., where Nis the number of equal space- 


intervals (or space-increments) into which the system is divided. 


The initial and boundary conditions are given as follows: 


Initial condition: p, = 0 @j =0 (6-14) 
J#i/2 jJt/2 d 
Inner boundary condition (predictor): Haste OW pe eee (6-15) 
2Ax d, 
Yall exdifiel| d 
Inner boundary condition (corrector): fas Loe (6-16) 
2Ax d, 
Outer boundary condition: p, > 0 as i > °%, for all (6-17) 
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The expressions for the predictor and corrector will now be written for the various grid 


points into which the system is divided. By defining 


ie (6-18) 


where Ax is the space increment, the predictor (Equation (6-12)) can be expressed in the 


following manner: 


n-1 


[24a fd, | | pee SR ile a A,(Ax)(d, / d,) 


n-1 


-A,(d, /d,)* op) -2Ax(d,/d,), fori=1 (6-19) 
n—1 
j+ul2 A; (i-1) Ax pi, - Di. ” j+il2 j+l2 _ (Ay)? A Da Day 
P3-) —|2+ Ae Toy tar a lpi + pi = (Ax) A, 2Ax 
es 
_Aedste| Pia — Pin)" opi, for 2 <i<N-1 (6-20) 
2Ax 
n—-I E ; 
s(N= jet — pi, : j+1/2 Pe fp 2 Prey — Dye 
pict 2+ atte Bhs Bls)* yp? pl? =(ax) (Pah Pl 
; ete 
Rrveaa | Pal PN, | lee ee tong (6-21) 
pAe: jad iE Lb Dy — Pi 1 f0rd 
A,é [ ae Pu — Pn+ 


Equations (6-19) through (6-21) represent a system of N equations in N unknowns 


(pressures), with the coefficients of the unknowns forming a tridiagonal matrix. The system 
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of equations can be solved by using the Thomas algorithm. The corrector (Equation (6-13)) 


may be written in the following manner 


-[2+(4 fd.) * | pit! + 2p! -|2-A4 /d,)* ap —2pi 
+2A,(Ax)"(d, /d,)—4Ax(d, /d,), fori = 1 (6-22) 
nol 
i+] _| 4 4 pAsli-Dax joy hdr — pi? \n Teles pita yl ee ele 
Di-} 2 eS ONDE Pigs = oy tee ~ Pis1 ) 
n-l 
(oi 5 2 A eds(i-Wax ae ee ae A j ete ray Vee) 
(D1, +pi.1) a A,e A Ax Cap, lone == (6-23) 
n-] 
j+l 2+A A;(N-1) Ax Dae a il = A Axl pitt? jrll2 j+l 
De tas 2€ oes O|\Dy =“, (pie — Pye )— pi, 
et j 2 — A eAs(N-Uax Dia SHA es i fori=N 
oe, 7 ine 2 A,e oe Q@ |Py, Or l = (6-24) 


Here also, Equations (6-22) through (6-24) form an NxN tridiagonal system of equations 


which can be solved by means of the Thomas algorithm. 
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Chapter VII 


ANALYTICAL SOLUTIONS FOR POWER-LAW FLOW THROUGH 
A FINITE-SIZED FRACTAL RESERVOIR 


In this chapter, analytical solutions of the partial differential equation (Equation (4- 
44)), that governs the transient flow of a non-Newtonian power-law fluid through a fractal 
reservoir, will be presented for a finite circular reservoir case. Both closed and constant- 
pressure outer boundary conditions will be considered. Moreover, the solutions will be 


presented for both constant-rate and constant-pressure inner boundary conditions. 
7.1 Constant-Rate Inner Boundary, Closed Outer Boundary 


The approximate (linearized) partial differential equation governing the radial flow of a 


non-Newtonian power-law fluid in a fractal reservoir is given by (Equation (4-36)) 


(d,/d, )(n+1)-d, (n-1)+n-3} 
On | jp CS RCA cae al ah ec 
Or? : TEE VON alltel Oath wit ote Ot 
with the initial condition given by 
p(r,0) = p,,forallr (7-2) 


For the constant-rate case, the inner boundary condition may be expressed as 


op 
or 


pe ads (7-3) 
k, | 2m, 


and for a closed (zero-flux outer boundary) reservoir, the outer boundary condition may be 


r=ly 


written as 
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op 
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where r, is the radius of the outer boundary of the circular reservoir. 


Before solving Equation (7-1), it is written in terms of dimensionless variables. 


Defining the dimensionless variables as follows 


(D.- Pp) k, te d, td, 


a een AS/ 50/571, eo) 
: 
i ee (7-6) 
|e 
and 
idea) tf 
D _ er (7-7) 


Equations (7-1) through (7-4) may be written, respectively, as 


Dp ‘ c hs ] d, a n Pp z (a, 1d,)° rXterisinev—artor-2} Po 


flee ere eat || 75 oF, Oty 
(7-8) 
Dp(iy.v) = O;for allrp (7-9) 
d 
CEO eri for all tr (7-10) 
Ce d, 


and 
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where 


212) 


Equation (7-8) may be solved, as in the infinite reservoir case, by applying the Laplace 


transformation. Application of the transformation to Equation (7-8) and the initial and 


boundary conditions yields 


aD, M d; d; aD p 2 {(d,/d,)(n+1)-d,(n-1) 
=F 4 Pe Sree ap = (d, /d,) r,| f ; IB, 


Ze 
dre i nd, d 
(7-13) 
d,/d 
TS) PE ee ea (7-14) 
ay |, l 
and 
Dy 25 (7-15) 
dry hig 


In order to solve Equation (7-13), it needs to be transformed into a Bessel equation by 
using a suitable change of variables (as has been shown in Chapter V). Following exactly 


the same procedure as outlined in Chapter V, it can be shown that the general solution to 


Equation (7-13) may be expressed as 
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Bein =e \C(ore) + C.K (ore )| (7-16) 


where 
Va) peers 
® n+1—d(n—1) ely) 
2V1(d, /d,) 
* (Glahart=d,n—)’ TaD 
72 d,in+J)—d,(n—I) 
deta I)-d,(n-1 
p= (eae ae rae 
and 
$f = Tg Ear aaa (J=20) 


2 


and where C, and Cy are constants to be determined from the two boundary conditions 


(Equations (7-14) and (7-15)). 


Differentiating Equation (7-16) with respect to rp and using the following recurrence 


relations (Carslaw and Jaeger, 1959), 

Ki(z) = -—K,(2)—K,_(2) (7-21) 
and 

I (2) = -"h(2)+1,4(2), (7-22) 


where z = arf, one gets 
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Applying the inner and outer boundary conditions to Equation (7-23), one obtains 


Ky_ eae 


~ PA{K,_(@l,(a78)— K,_.(ar8)h ala] aor 


and 


Leh (ars, ) 
P?1K,_,(0t)L,-a( Oren) - K,_.( arp), (o0)] 


(7-25) 


2 


Substituting for C; and C> in Equation (7-16), the pressure solution is obtained as 


ri Bil (ors, )K, (cr) +1, (crf )K,_,(ar'5)} 


We | 
Pool! = “pal, (al, (ard) -K, ,(or8), 0} (7-26) 


where v, @, Band yare given by Equations (7-17), (7-18), (7-19) and (7-20), respectively. 


Equation (7-26) is the general solution (in Laplace space) for the transient pressure 


behaviour in a circular closed fractal reservoir with a centrally located well producing a 


power-law fluid at a constant rate. Equation (7-26) reduces to the following form at the 


wellbore (7p = 1) 


1,_,(c8, \K, (a) +1,(@)K,_,(7r5 ) 


+. 
; _ __ Falaro)K(@)+(@)K (or) 7-27 
Pyp(l) prtK POON (are )- Kuss Orne (a)} ( ) 
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For d¢= d, = 2, Equation (7-27) exhibits the wellbore pressure response during the flow of a 


power-law fluid through a closed (circular) homogeneous reservoir (Equation (6-16) of 


Ikoku (1978)). 
7.2 Constant-Rate Inner Boundary, Constant-Pressure Outer Boundary 


This case represents the situation where radial flow takes place in a circular fractal 
reservoir, the outer boundary of which is essentially at a constant pressure due to either 
natural or artificial pressure maintenance at that location. Here also, as in the previous case, 
the objective is to obtain the transient pressure solution of Equation (7-1) with the initial and 
inner boundary conditions being expressed by Equations (7-2) and (7-3), respectively. In 
this case, however, the outer boundary condition is different from that in the previous case 


and may be written as 
Dita tio—p; (7-28) 


In a manner similar to that used in the previous case, one can rewrite the governing 
equation (Equation (7-1)) in dimensionless form (where the dimensionless variables are 
defined by Equations (7-5) through (7-7)) and then express the resulting equation and the 
boundary conditions in Laplace space. The Laplace transform of the governing equation in 
dimensionless form is given by Equation (7-13); the Laplace transform of the inner 
boundary condition is given by Equation (7-14). The Laplace transform of the 


dimensionless form of the outer boundary condition is expressed as 
Pp("p.1) = 0 (7-29) 
Equation (7-13) has the following general solution 


Dy(tp,l) = rE [Cy,(crp) + CK,(075 )) (7-30) 
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where the parameters v, a, B and yare given by Equations (7-17), (7-18), (7-19) and (7- 
20), respectively, and C, and C> are constants to be determined from the boundary 


conditions. Differentiating Equation (7-30) and using the recurrence relations given by 


Equations (7-21) and (7-22), one obtains 


oo = Creal GLa Ors )-OK,_,(arp)) (7-31) 


Using the inner and outer boundary conditions, given by Equations (7-14) and (7-29), 
respectively, the constants C, and Cy can be determined; substituting for C, and Cy in 


Equation (7-30), the pressure solution can be expressed as 


P 
Pp(%,l) = ra (7-32) 


Equation (7-32) is the Laplace transform of the general solution for the transient pressure 
behaviour for the case of constant rate of flow of a power-law fluid in a circular fractal 


reservoir with a constant pressure outer boundary. At the wellbore (rp = 1), Equation (7-32) 


reduces to the following form 


B - B 
Dap (1) = aoa I,(a)K, (ard ) (7-33) 


K,_,(oe)1, (or8,) + K, (ares )L,.(@)} 


For de = d,= 2, Equation (7-33) reduces to Equation (6-28) of Ikoku (1978) for a 


homogeneous reservoir. 


7.3 Constant-Pressure Inner Boundary, Closed Outer Boundary 
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The linearized partial differential equation governing the transient flow of a power-law 
fluid in a fractal reservoir is expressed by Equation (7-1). Before solving Equation (7-1) for 
a variety of boundary conditions, it is expressed in terms of dimensionless variables. For the 


constant-pressure inner boundary condition, the dimensionless pressure may be defined as 


follows 


n= JON Iaith 
Dol Tn tpd = Bi - P(r,t) (7-34) 
P; ee bee 


where py, is the constant wellbore pressure. The dimensionless radius and time are 
expressed, as in the previous cases, by Equations (7-6) and (7-7), respectively. It can be 
shown that, by using these definitions of dimensionelsss variables, Equation (7-1) can be 
expressed in dimensionless form as Equation (7-8). Assuming that initially the pressure 
throughout the reservoir is uniform (such that @ tp = 0, pp = Oat all rp), Equation (7-8) 
can be rewritten, using the Laplace transform, as Equation (7-13), which can be transformed 
into a modified Bessel equation and solved in a straightforward manner. The general 


solution to Equation (7-13), as has been shown earlier, is of the form 


Dorrie .|Clh(a7, ) ot C,K,(ors }| (7-35) 


where the parameters v, a, B and yare given by Equations (7-17), (7-18), (7-19) and (7- 


20), respectively, and C; and C> are constants to be determined from the inner and outer 


boundary conditions. 


For a flow situation corresponding to this case, the inner and outer boundary 


conditions are expressed, respectively, as 


p(r=r,,t) = PD, (7-36) 
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and 


op 
rr () : 
Ble, (7-37) 
where r, is the radius of the outer boundary of the circular reservoir. Rewriting Equations 


(7-36) and (7-37) in terms of dimensionless variables, and then using the Laplace transform, 


the inner and outer boundary conditions are expressed as 


P- Ji 

Pp(% =1,1) = 5 (7-38) 
and 

Dp a (7-39) 

dr, set 


where rp is the dimensionless radius of the outer boundary of the reservoir (given by 


Equation (7-12)). 


Using Equations (7-35), (7-38) and (7-39), the constants C, and C2 can be determined 


as follows 
B 
C= K (a) ; (7-40) 
NK, (co), .,(o75,) + K,_,(ar5)L,(@)} 
B 
one foe rn) (7-41) 


and thus, the pressure solution becomes 
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Pp(%,!) = (7-42) 
NK, (orl, (ar) + K,_,(ord\1,(ce)} 
Defining the transient dimensionless rate as 
ee are Peg (7-43) 
°K, (B=) dy td, \ ahr, )’ 
where the dimensionless rate is related to the pressure gradient at the wellbore as 
qp(d,/d,) = - Pp (7-44) 
'D Tp=l 


and taking the Laplace transforms of Equation (7-44) and using Equation (7-42), one 


obtains the following expression for the dimensionless rate in Laplace space 


I,_,(ar', )K,_,(@) a Ke (orn \I,(@) 


a (7-45) 
ATK, (ODE, (78) + K,_, (0755 )L,(@)} 

Defining the cumulative production (over a given time fp) as 
Gy = |. ap diy (7-46) 


and taking its Laplace transform, the cumulative production solution in Laplace space can be 


expressed as 


¥ fe oor \K _,(a)—K,_,( cr I, (0) 
PAK (cel, (075) + K-( Oren) ()} 


(7-47) 
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7.4 Constant-Pressure Inner Boundary, Constant-Pressure Outer Boundary 


In this case, the inner and outer boundary conditions are defined by Equations (7-38) 
and (7-29), respectively, with the dimensionless pressure being defined by Equation (7-34). 
Using the two boundary conditions, the two unknown constants in Equation (7-35) can be 


determined and the dimensionless pressure solution expressed as 


age cers \K, (org) — K,( arg I, (org, )} 


Plo ~~ Harb \iCa) Kai (or8)} 


(7-48) 
and thus, the dimensionless rate and cumulative production solutions for this case are 


I,_,(@)K, (ork, ) + K,_,(a)I,(arg) 


00) = “TilK,(c, (orb) —K,(ar8,)1,(co} 


(7-49) 
and 


ee (a)K, (orr', )+ K,_,(@), (ars, 


Q,(1) = 74K, (oe)1, (or) = K, (orf, \,(a)} 3 


(7-50) 


respectively, where the dimensionless rate and cumulative production are defined by 


Equations (7-43) and (7-46), respectively. 
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Chapter VII 


ANALYTICAL SOLUTIONS FOR POWER-LAW FLOW THROUGH 
A TWO-ZONE COMPOSITE RESERVOIR 


In the realm of hydrocarbon reservoir engineering and well testing, the study of 
composite reservoirs is of great importance because of the wide range of reservoir 
configurations they represent. The composite reservoir model has been used to describe and 
study the behaviour of damaged/stimulated reservoirs, reservoirs undergoing waterflood, in- 
situ combustion, and so forth. In a (two-zone) composite reservoir model, the reservoir 
system is generally considered to include a circular inner region with rock and fluid 
properties significantly different from those in the outer region. It is usual in these modelling 
efforts to assume that the radial extent of the inner zone and the flow properties (such as 
permeability) of the two zones are adequately described by the parameters of the composite 
model. 

In this section, an infinite reservoir with a physical radial discontinuity in the rock 
system is considered. The zone nearest to the wellbore is a fractal reservoir (with spatial 
distributions of porosity and permeability); the outer zone, infinite in extent, is 
homogeneous. The major contribution of the material presented in this section is the 
development of analytical transient pressure and rate solutions during the flow of a non- 
Newtonian power-law fluid in such a reservoir system. The characteristics of the solutions, 
presented in the form of dimensionless pressure and pressure-derivative curves, will be 


discussed in a subsequent chapter. 
8.1 Mathematical Formulation 


The reservoir system under consideration is illustrated schematically in Figure 8-1, 
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Figure 8-1: Schematic Representation of a Two-Zone 
Composite Reservoir 
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which exhibits a reservoir system composed of two concentric zones. The inner zone is 
composed of a fractal zone; that is, the inner zone consists of a fractal network of fractures 
embedded into a homogeneous matrix with the former dominating the flow, or, 
equivalently, the entire inner region is fractal. The inner zone surrounds the well to a radius 
r, and there is a sharp radial discontinuity between the inner and outer zones. The outer zone 


is homogeneous, has a permeability k’ and a porosity @’ and is infinite in extent. The 


wellbore of radius ry, is located at the centre of the reservoir system. 


The transient flow behaviour of a single-phase, slightly compressible, power-law fluid 
in the reservoir system described above and shown in Figure 8-1 is considered. The 
transient cylindrical flow in the composite system is modelled by the following system of 


partial differential equations: 


Inner zone 


(d,/d,)(nt+1)—d, (n—1)+n-3 
Da . E ermrates ee a BO 
f 


Or? Heatoweed. |r or Tae ie Ot” 
LOL a (8-1) 
Outer zone 
oP + a = oo forr > rj] (8-2) 
where G is given by Equation (4-37). Also, 
ae I-n n=l 
G’ = ae [Zossinrasoe yp (4) (8-3) 
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and p, and p2 are the pressures in the inner and outer zones, respectively. Transforming 


Equations (8-1) and (8-2) into dimensionless form yields 


2 
CLIT li eee Wc a (dy, ) rphldrt@ert-arin--2} Pos 
Ory eek ter Aes, Cole Seg ee Oty 
foray a (8-4) 
rape n OP p> 1-n OP 
=a ae aa = (d, /d,) Onl. a forrp >a (8-5) 
where 
re tN 8-6 
: (8-6) 
and 
l=-n 
Lae bee oe 
2, = AE tal (8-7) 
9,6, 1k, \ ,K,, 
The dimensionless variables are defined as follows 
(Veh )k, as d / d, 
ge oem EE a (8-8) 
A (q/2mh) 
ae 16 (8-9) 
fs 
2 
i (8-10) 
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8.1.1 Constant Rate Case 


Equations (8-4) and (8-5) are first solved for a constant rate inner boundary condition 
and an infinite outer boundary condition. The various boundary conditions for the flow 


situation considered in this case are described below. 
Inner boundary condition: 


This condition is for a finite-sized wellbore producing at a constant rate such that 


Por 


Bi for all t (8-11 
— =, or a 
Ory : d D ) 


=] Ss 
Interface boundary conditions: 


Two interface boundary conditions will be considered. The first condition is used to impose 


pressure continuity at the interface between the two zones: 
Pp(4,tp) = Pp2(A,tp) (8-12) 


The second condition imposes rate continuity at the interface such that the flow rate from the 


outer zone across the interface equals the flow rate into the inner zone: 


Por = EPA (OGRE (8-13) 
Or, Ig 
where 
l-n 
og ea aa (8-14) 
k,\ o’k 
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Outer boundary condition: 
For an infinite outer boundary condition, one has 
Pp2 > 0 as tH > (8-15) 


The system of partial differential equations and the associated boundary conditions, as 


shown above, is solved using the Laplace transformation. 


Applying the Laplace transformation to Equations (8-4) and (8-5) yields 


; a’Dp, bs "| 444), Gu = (ayia lrg eae alone 
n 'D 


(8-16) 


25 2 
pF Por 5 py, Por - (a sidy\copramhps: (8-17) 
ary lp 


It is to be noted that in transforming Equations (8-4) and (8-5) into Equations (8-16) and (8- 


17), respectively, use has been made of the initial condition (i.e., @ tp = 0, pp = 0). The 


boundary conditions become, in Laplace space, 


Dor Ste aie (inner boundary condition) (8-18) 
dr, es l 

Dp, (4,1) = Dp2(a,l) (pressure continuity at interface) (8-19) 
QD = o, Poa @ rp = a (rate continuity at interface) (8-20) 
dr, ies 
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and 
Pp2 7 9 as Hy (8-21) 
As has been demonstrated in Chapter V, the solution of Equation (8-16) is of the form 


Boi(tp.l) = rp.[C,I, (arg) + C,K,(arp )| (8-22) 


where v, a, Band yare given by Equations (5-21), (5-17), (5-18) and (5-19), respectively. 


In a similar fashion, it can be assumed that the solution of Equation (8-17) is of the form 
Poo = "5 8(p') (8-23) 


where 


, 


Wee (8-24) 


Following an approach identical to the one outlined in Chapter V, it is possible to obtain, by 


substituting Equations (8-23) and (8-24) into Equation (8-17), the following equation in p' 


fee 
pre” + p's’ -|5—"|2- p%e = 0 sige? 


In arriving at Equation (8-25) it was necessary to choose 


if eee he (d, /d,) (8-26) 
3-n 
3-n 
ine 8-27 
B 5 (8-27) 


and 
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Equation (8-25) is in the form of the modified Bessel equation and, hence, the solution to 


Equation (8-17) is 


Pp2(tp.!) = rh |Csl,-(ar rf ) + CyKy (crf ) (8-29) 

where 
 apleiate 8-30 
par ( si ) 


By applying the outer boundary condition (Equation (8-21)) to Equation (8-29), it is easily 


seen that C3 = 0 and thus, Equation (8-29) reduces to 
Bo2(%.l) = C, 1p K,(0 rp ) (8-31) 
From Equations (8-18) and (8-22), the constants C, and C> can be shown to be related as 
Gn OC pC. Kw (a= 1s” (8-32) 


Substituting Equations (8-22) and (8-31) into the pressure continuity relationship (Equation 


(8-19)) yields 
Gon C.K (on) = Gra" K,(o'a" ) (8-33) 


Finally, by applying Equations (8-22) and (8-31) in the rate continuity relationship 
(Equation (8-20)), one obtains 
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Solving Equations (8-32) through (8-34) simultaneously, the three unknowns C, , C> and 


C4 can be determined. The constants C and C> are found to be 


AK,.(a a” )K,_,(aa®) — K,(aa®)K,_.(a'a®’) 


(oe — A (8-35) 
and 
Fyn 5S” B B off 
C, = Ahy(ola He (aa") + 1,(oa" )Ky (oa) (8-36) 


A 


respectively, where 


0,0, (8-37) 
and 
A = PAK, (aa® {1,_,(oa")K,_,(a) — 1,_,(00)K,_, (aca? jh 
+K,,(o'a” {I,(aa")K,_(0) + 1,..(00)K,(aa" )}] (8-38) 
Thus, the required pressure solution at the well in the inner zone is given by 
Pp (1,1) = Cl,(a) + C,K,(a) (8-39) 


where Cj and C> are given by Equations (8-35) and (8-36), respectively. One may write the 


above equation also in the following fashion 


Dp (Ll) = a (8-40) 
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where 


a, = {AK,.(o'a® )K,_(aa") -— K,(aa")K,_,(o a® by (cx), (8-41) 


a, = {AK,(ara® )I,_(aa’) + 1,(ca® )K,_,.(or a” )\K,(a) (8-42) 


and A is given by Equation (8-38). 
8.1.2 Constant Pressure Case 


It has been demonstrated in the previous section, by using the initial and the outer 
boundary conditions, that the pressure solutions for the inner and outer zones are given, 
respectively, by Equations (8-22) and (8-31). For the constant pressure inner boundary 
condition, defining the dimensionless pressure as 

P= P 
[ikke SS aes (8-43) 
PD; = Py 
and using the same initial and outer boundary (infinite) conditions as in the previous case, 
the pressure solutions for the inner and outer zones may again be expressed by Equations 


(8-22) and (8-31), respectively. 


The inner boundary condition in Laplace space is 


= I : 
Poily= ~ 7 (8-44) 
Substituting Equation (8-22) into Equation (8-44), one gets 
: 8-45 
Gao C,K.(a) = 7 (8-45) 
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Applying Equations (8-22) and (8-31) in the pressure and rate continuity relationships, one 
obtains, as in the previous case, Equations (8-33) and (8-34). Solving these two equations 
together with Equation (8-45) simultaneously, it is possible to determine C}, Cy and C4. The 


constants C, and C> are determined, respectively, to be 


AK,.(at' a” )K,_,(aa®) - K,(aa®)K,_,(a’a® ) 


C, = 
L{a,+ a} 


(8-46) 


oe AK, (a'a® )I,_ aa’) + I,(aa’ )K,_,.(a’a* ) (8-47 
: Via, ai) 1) 


where a7 and ap are given by Equations (8-41) and (8-42), respectively. Thus, the pressure 


solution for the inner zone is 
Poiltp.l) = rB[Cyly(arp) + CoK,(arp)| (8-48) 
where C and C are given by Equations (8-46) and (8-47), respectively. 


Defining dimensionless rate and cumulative production, respectively, as 


Ar q é 
= ___ 8-49 
ve Ke (D; — p,,) d; / d, (2 


and 


Oy SCS, (8-50) 


the dimensionless rate and cumulative production can be evaluated in Laplace space as 
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Gp(d, /d,) = pa (8-51) 
D \Wrp=l 
and 
One= 4 (8-52) 


The expressions for the dimensionless rate and cumulative production can thus be written as 


=) es ae 4 
Gp(l) = P (a, +4,) (8-53) 
and 
= A 
St eal aa -54 
0, (1) ie (a,+a,)’ (8 5 ) 


respectively, where A, a; and a are given by Equations (8-38), (8-41) and (8-42), 


respectively. 
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Chapter IX 


PRESSURE TRANSIENTS WITH MATRIX PARTICIPATION 
IN FLOW: A SPECIAL CASE 


The previous chapters dealt with the development of transient pressure and rate 
solutions for the single-phase flow of a non-Newtonian power-law fluid through infinite and 
finite fractal flow networks. In the material considered in the previous chapters, the analysis 
of the pressure transients was performed by assuming that the matrix does not participate in 
the flow process. In this chapter, a mathematical formulation is presented for the transient 
flow of a non-Newtonian power-law fluid in a naturally fractured reservoir that consists of a 
fractal network of fractures embedded into a homogeneous matrix. The conventional 
Newtonian flow model in a double-porosity system is generalized by allowing the matrix to 
exchange fluid with the fractal fracture system. Following the Warren and Root approach 
(Warren and Root, 1963), it is assumed in the present model that flow between the matrix 
blocks takes place only through the fracture system. Furthermore, the well intersects the 
fracture system and all the fluid produced at the wellbore moves through the fractures. 
Analytical solutions are obtained for a special case of the flow situation where a Newtonian 
fluid travels through infinite and finite fracture/matrix systems; the main reasons for 


considering this particular case are also discussed. 
9.1 Mathematical Formulation 


The continuity equation for cylindrical flow in this fracture-matrix system may be 


written as 


ed 


a) 
(rp,u,) = <(6,p,) 2 5 (Onin) (9-1) 


a 
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where the subscripts f and m refer to the fracture network and matrix, respectively. The 


modified Darcy's law for a power-law fluid flowing in the fracture network is 


si k,(r) op, In 
u, = [a ir) | (9-2) 


Applying Equation (9-2) in Equation (9-1) yields 


Ops By il en ene es 
k Sig fy ead een a 
ma ae ie el a 2 Nibg(r) dr | 


Op OD,, l/n op 
= |,(rde, Be ma eeG,. nt Bauer Heir e| (9-3) 


In deriving the above equation, it was assumed that pressure gradients in the fracture 
network are small at all times. Equation (9-3) can be further simplified, by making use of 


Equations (4-13), (4-14), (4-16) and (4-34), to 


or? uy n nd. a. |r x 
aa o,c,(rir,)! aaa PE it! 
ert or or le (9-4) 


For a Newtonian fluid (n = 1) flowing in a naturally fractured reservoir where a 


homogeneous fracture network is embedded into a homogeneous matrix (such that dr = d, 


= 2), Equation (9-4) is reduced to the following form 


of oh tye 


OP, , 15, _ a wae a (9-5) 
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Equation (9-5) is the Warren and Root (1963) equation for transient radial flow in a 
double-porosity system (Sabet, 1991). Warren and Root have shown that, for such a flow 
medium, a semilog plot of the drawdown data (flowing wellbore pressure versus time) 
exhibits two parallel straight lines. The earlier line indicates transient radial flow through the 
fractures before the matrix makes its presence felt; the later straight line develops after an 
equilibrium is reached between the fracture and the matrix pressures. The transition between 
these two lines develops as a result of the matrix-to-fracture interporosity crossflow. This 
transition from early production from the fractures to late production from the total reservoir 
(matrix and fractures) is affected by the way the matrix and fracture network are assumed to 
interact. In the Warren and Root model (Warren and Root, 1963), the flow from matrix to 
fractures is assumed to take place under pseudosteady-state conditions; in other words, the 
interporosity flow rate is proportional to the pressure difference between matrix and 
fractures. Such an assumption was also employed by Chang and Yortsos (1990) in their 
model describing flow in both the fractal object and the matrix. The pseudosteady-state flux 


assumption will also be used in the present study. 


Following the approach of Chang and Yortsos (1990), the expression for the 


interporosity fluid exchange rate is given by 


l/n 
Gas ie lB) (9-6) 


ewer 


where D’ is the fractal exponent for the perimeter of the fractal object and b is the 
corresponding proportionality constant with a dimension of [L?~?], and D” is the fractal 
exponent for the average distance between the matrix and the fractal object and e is the 
corresponding proportionality constant with a dimension of [L!~?’]. In the Euclidean limit, 


D’'=2,D" = 0, b = 1 ande = the characteristic length / of the Warren and Root model 
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(Chang and Yortsos, 1990). By applying a mass balance for the fluid contained within the 


matrix, it can be shown that 


] 


ff — Pm me 
= (PnP mn) (9-7) 


2mrh 


Combining Equations (9-6) and (9-7) and simplifying, one obtains 


l/n 


ap,. pr D’-2-D"In 12k,,( D, S. Dn) 


ote" 2thd,.c 


™ | H(9+3/n)"(150k,o,) 2 


Equations (9-4) and (9-8) can be combined and solved for various sets of boundary 
conditions to present and analyze the behaviour of the pressure transients with matrix 


participation in flow. 
9.2 Analytical Solutions for a Special Case 


Equation (9-4) is an approximate partial differential equation describing the transient 
flow of a non-Newtonian power-law fluid through a fractal fracture network-matrix system. 
The suitability of Equation (9-4) in describing such flows is limited mainly because of the 
use of the approximation, given by Equation (4-29), in deriving it. However, for the 
Newtonian fluid case (n = 1), Equation (9-4) does not suffer from this limitation. It is, 
therefore, of theoretical and practical value to consider analytical solutions for Equation (9-4) 
for the case of a Newtonian fluid. For the Newtonian fluid case, Equations (9-4) and (9-8) 


reduce, respectively, to 
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and 


OD n bk,, 
ot eu 2Thd,,.C, 


-O 


(p,-D,)r (9-10) 

where 0 = D” + 2 —D’. The similarity between the above equations and the corresponding 
equations (for flow in both fractal object and matrix) used in the study of Chang and 
Yortsos (1990) is obvious; in fact, Equation (9-10) is identical to Equation (25) of Chang 
and Yortsos for d = 2. In their study, however, Chang and Yortsos considered numerical 
solutions to the partial differential equations governing flow in both fractures and matrix; 
analytical solutions could not be considered because of the presence of spatially variable 
coefficients. It was found in their study that the system pressure response exhibited early- 
and late-time linear behaviour in the log-log plot (of dimensionless pressure against 
dimensionless time). It can be inferred from the results presented in their study that the slope 
of the early linear segment depends on d,, whereas that of the late segment depends on both 
drand d,; it was, however, mentioned that the difference in the two slopes is small. Chang 
and Yortsos also studied the effects of the various relevant parameters on the transitional 


period between the two linear segments and concluded that ohad the least significant effect. 
9.2.1 Constant Rate Case 


Defining dimensionless pressure, time and radius as follows 


_ 20k, h (dy! d,) (Di= Pym) (9-11) 
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T= were (9-12) 
i (9-13) 
where 
pee ro. (9-14) 
Pyle + OmCm 


It is to be noted that the dimensionless time defined here (Equation (9-12)) and that defined 


in the previous chapters (Equation (4-40)), for n = 1, are related as T= Wtp. Equations (9- 


9) and (9-10) can be expressed in terms of the dimensionless variables as 


2 2 
0 Pog zm 4, 1-7) 1 Su = (= jas OP py aye) OP om oe 
d 
D s 


Cig i UG d, OT oe 
(9-15) 
and 
Eee AA (py | Eup) (9-16) 
OT igs 
where 
Dp=D" 
bk, (9-17) 


~ (1— @)2mhek, (d, /d,)? 


The parameter @ represents the ratio of near-wellbore fracture storage to total storage. The 


interporosity interaction parameter, A, is proportional to the ratio of matrix permeability to 
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near-wellbore fracture permeability (k,,./k,,) and generally has values much smaller than one. 
Large values of @ (say, greater than 0.01) would indicate a significant amount of storage in 


the fracture system; large values of A (say, larger than 10-5) would indicate a relatively large 


flow conductivity in the matrix. 


Before solving the equations describing the transient response of the fracture-matrix 
system, two assumptions will be made — the system behaviour for dy = 2, and for o= 0 will 
be considered. The effect of the second assumption will be felt only during the transitional 
period, as has been discussed earlier. This effect, however, will be minor, as the parameters 
@ and J have been shown to have more significant roles to play during that period (Chang 
and Yortsos, 1990). The first assumption is extremely significant; it will influence both the 
transitional period and the late-time behaviour of the system. However, as will be shown in 
the next chapter, this assumption is intended to help bring out some interesting 


characteristics of the transient response of the fractured medium. 


Incorporating the above-described assumptions, Equations (9-15) and (9-16) can be 


expressed in Laplace space and then combined to yield 


da’p Ae dp oe 
ea te 9-18) 
D S 
where 
E = 1(2/d,) (a4) (9-19) 


The inner and outer (infinite) boundary conditions are 
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Ppp 70 as > (9-21) 


It is assumed that the pressure solution for the given drawdown problem may be expressed 


as 
She SIRI ene (9-22) 


Substitution of Equation (9-22) in Equation (9-18) yields 


per’ + pf [1-2] f- pif = 0 (9-23) 


provided one chooses 


rate 7 i ware (9-24) 
aa 


B= 2/d, (9-25) 
and 

y = 2/d, - 1 (9-26) 
The solution to Equation (9-23) is, as has been shown previously in Chapter IV, 

f(p) = Cl,(p) + ©,K,(p) (9-27) 
where 


y=1-d,2 (9-28) 
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Combining Equations (9-22) and (9-27) and applying the outer boundary condition 
(Equation (9-21)), one gets 


Po, (Mp1) = Cy 0h K, (org ) (9-29) 


Applying the inner boundary condition (Equation (9-20)) in Equation (9-28), the constant 


C is determined to be 


Cr ea eal 9-30 
ak (a) eu) 


and thus, the dimensionless wellbore pressure solution in Laplace space is obtained as 


Ky (a) 


el) er 
l-v 


(9-31) 
For the case of a closed circular reservoir with a centrally-located well producing at a 
constant rate, the wellbore pressure solution may be obtained in a similar fashion 


I,_,(0re )K,(@) +1,(0)K,_,(cr) 


y-l 


Pie Rca on coer ead oe (9-32) 
al{K, ,(@)L,.(arg)-K,_,(a75 )L,(@)} 
The wellbore pressure solution for a constant pressure outer boundary condition is 
1,( cr’, \K,() - 1,(a)K,( ort 
Pp, (% =1, 1) = | 2) | 2 (9-33) 


od} K,_,(a)I,( 75 ) +K, (ard, )I,_.(a)} 


It can be shown from Equation (9-24) that as a1, A> vi and Equations (9-31) 
through (9-33) reduce to the corresponding wellbore pressure solutions for flow only in the 


fractal fracture network for n = 1 and d= 2. This behaviour is expected, since the reservoir 
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must contain only the fracture system if @ is to approach unity. Similarly, Equations (9-31) 
through (9-33) indicate fracture flow behaviour as 2 > oo; this is also proper because either 
@— 1 or k,, - © (indicating no impedance to interporosity flow) if A has to approach 


infinity. 
Limiting Solutions for Early and Late Times 
At early times, 1 4 0, so that @ > Val and the modified Bessel function, K,(z), 


may be approximated by Equation (5-31). Thus, Equation (9-31) approximates to 


K, (Val) 


Po (% = 1,1) = Tal K,,. (lal) (9-34) 


Using Equation (5-31), the early-time limiting wellbore pressure solution in the real plane 


can then be expressed as 


Pog (% = ], T) = 24) = 2,|2 (9-35) 


where Tis given by Equation (9-12) and tp by Equation (4-35), for n = 1 and dy= 2. Thus, 
at early times, the wellbore pressure response would be exactly the same as that when only 


the fracture system participates in the flow process (see Section 5.1.2). 


At late times, ] > 0 (a > VI), and applying Equation (5-50) in Equation (9-31), the 


late-time approximation of the Laplace space wellbore pressure solution is 


oo iy) 


sr (9-36) 
ITd-v)a@ 


Pog (tp = 7 1) = 


which can be inverted to 
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(9-37) 

It is also of interest to consider the early- and late-time limiting solutions for the special case 
of @ = 0, which applies for negligible storage capacity in the fracture system; this case has 
been discussed extensively for a homogeneous fracture system by Warren and Root (1963). 
For this case, at early times, when a > V1 , the wellbore pressure solution has a constant 


value given by 


K, (va 
Pos (% =a pe (9-38) 


Thus, a constant early-time wellbore pressure in a "double-porosity" system would indicate 
an extremely small value of the storage parameter, @. At large times, @ > vi, and the 


wellbore pressure response is given by Equation (9-37). 
9.2.2 Constant Pressure Case 


For the constant pressure inner boundary condition, the dimensionless pressure has 


been defined previously (Equation (5-77)) as 


Se f BM 
Lf ma a 


P= Py 


(9-39) 


The dimensionless time and radius are given by Equations (9-12) and (9-13), respectively. 
Equation (9-18) describes the transient flow of a Newtonian fluid in a fractal fracture- 
homogeneous matrix system with dr= 2 and o= 0. The inner boundary condition is 

t 


Hil, Valet ee (9-40) 
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and the outer boundary condition (for an infinite system) is given by Equation (9-21). Ina 
manner similar to that shown for the constant rate case, the dimensionless fracture pressure 


solution for the constant pressure inner and infinite outer boundary conditions may be 


obtained as 
= eens) 
P(r, 1) = > : 
p(T» !) 1K, (a) (9-41) 
Defining the dimensionless rate as 
quid, 
— 9-42 
10 Gak,WD;—P,) ae 
it is possible to show that 
d OP ny 
<j =) Gated a 9-43 
Ip Sa (9-43) 


rp=l 


Combining Equation (9-43) with Equation (9-41), the following expression for 


dimensionless rate is obtained 


Gp a OR, (a) (9-44) 
IK, (0) 
Finally, by defining the cumulative production as 
05 = [andr, (9-45) 


it can be demonstrated very simply, by combining Equations (9-44) and (9-45), that 
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Q, = PR (a) (9-46) 


Similarly, the dimensionless rate and cumulative production solutions for a closed 


reservoir are obtained as 


onl, (on K,_,(@)-K,_,(ard)h,_.(01)} 


7-(l) = : 
De NK, (@)l, (ar) +K,_,(orf )1,(a)} Rage 
and 
aie a s(ors (2) - K(at M,.(ce)} ron 
1K, (aL, s(ar) + K,_.(are ),(a)} 
Finally, the solutions for a constant pressure outer boundary situation are given by 
P ofl, (0K, (org )+ K,_,(@)L, (on )} 
l) = — 9-49 
010 * ~ iT (ai, (ar) — (ax Jia} We 
and 
ne oe{1,_, (CLK, (or )+ K,_,(0e)1, (Orr )} non 


P{K,(a)I,( orf) K (ores t,o) 
Limiting Solutions for Early and Late Times 


At early times, / > o, so that & > 4/@l and thus, from Equation (9-44), 
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Equation (9-51) can be inverted analytically to yield 


@ 1 
Tp Oe EY ares Ey ar 
UT Mp 
and thus, the cumulative production solution becomes 


Op(tp) = 2,,-2 
ING 


At late times, 1 0 (a > V1), and thus, from Equation (9-44), 


page) 
D T(v) a2v-1 pi-v 


The late-time dimensionless rate solution, expressed in the real plane, is then 


q1-2v 
t = 
qp\tp) o’ Tv) tp 


and the cumulative production solution is 


(air 


Op(tp) = 2-1 Tw) (1) 
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(9-51) 


(9-52) 


(9-53) 


(9-54) 


(9-55) 


(9-56) 
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Chapter X 


RESULTS AND DISCUSSION 


Various models have been presented in the previous chapters to analyze the transient 
pressure and rate behaviour of power-law fluid flow in reservoirs exhibiting power-law 
variations of permeability and porosity with distance from the wellbore. In Chapter IV, a 
partial differential equation was derived for an approximate description of the transient flow 
of a non-Newtonian power-law fluid in an infinite fractal flow medium. Chapter V presents 
analytical solutions, in Laplace space, to this equation for both constant-rate and -pressure 
inner boundary conditions; early- and large-time limiting forms of the analytical solutions 
have been derived in real space. In Chapter VI, a finite-difference scheme was presented to 
solve the nonlinear partial differential equation describing the transient flow of a power-law 
fluid in an infinite fractal reservoir. Chapter VII contains analytical pressure and rate 
solutions for finite reservoirs; both closed and constant-pressure outer boundary conditions 
have been considered. Chapter VIII presents analytical pressure and rate solutions for a two- 
zone composite reservoir, with the inner zone being fractal and the outer homogeneous. 
Chapter IX demonstrates the development of analytical solutions for a special case of the 
situation where the transient flow of a Newtonian fluid takes place in a fracture/matrix 
system, with the fracture network showing anomalous diffusion. 

This chapter presents the results obtained from the different models described in the 
previous chapters (see Figure 10-1) and discusses the effects of the various parameters on 
these results. The discussion mainly focusses on analysis of results obtained from the 
transient pressure and rate models for flow in an infinite fractal medium, with particular 
emphasis on the late-time behaviour of the flow system. Discussion is also presented on the 


results obtained for flow in a finite system, in a composite system and in an infinite fractal 
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TRANSIENT FLOW IN A FRACTAL MEDIUM 
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Figure 10-1: An Organizational Map Depicting the 
Situations Considered in this Study 
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fracture/matrix system. 
10.1 Transient Flow Through a Fractal Medium 
10.1.1 Analytical Solutions: Infinite Reservoir 


Transient solutions for dimensionless pressure (Equation (5-30)) and rate (Equation 
(5-82)) for an infinite fractal reservoir can be evaluated readily in real space by using 
numerical techniques. In this study, the numerical Laplace transform inversion scheme 
developed by Stehfest (1970) is applied to obtain the dimensionless pressure (pyyp) and 
pressure-derivative (dpyp/dintp) solutions for various combinations of ds and n. Some of 
the results for d; = 2 and for various values of the flow behaviour index, n, are shown in 
Figures 10-2 through 10-4. It can be seen from these figures that with increasing time, the 
pressure solutions start diverging from each other depending on the relative magnitudes of 
n; the dimensionless pressure at any given time is larger for a smaller value of n. In fact, 
Figure 10-2 shows that the dimensionless pressure at a given time during the production (or 
injection) of a pseudoplastic fluid may be more than an order of magnitude larger than that 
for a Newtonian fluid. However, it is to be noted from Equation (4-38) that the definition of 
dimensionless pressure involves flow rate and well radius raised to powers of n (where 0 < 
n < 1) and thus, the actual pressure drop (for production) or increase (for injection) may be 


lower for a pseudoplastic fluid than for a Newtonian one (Ikoku, 1978; Olarewaju, 1992). 


Figure 10-2 shows that at late times and for v > 0, the dimensionless pressure solution 
exhibits a straight line in the log-log plot; this is expected, as has been demonstrated earlier 
by Equation (5-47). The slope of this straight line is equal to the parameter v. It is interesting 
to note, however, that the pressure-derivative plots (see Figure 10-3) exhibit a large-time 
straight line behaviour for a longer period of time; such a behaviour may be explained by 


inspecting Equations (5-43) and (5-47). The slope of the linear segment of the pressure- 
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Figure 10-2: Dimensionless Wellbore Pressure versus 
Time for an Infinite System; Pseudoplastic Fluids, d,s = 2 
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Figure 10-3: Dimensionless Wellbore Pressure 
Derivative versus Time for an Infinite System; 
Pseudoplastic Fluids, dg = 2 
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Figure 10-4: Dimensionless Wellbore Pressure Derivative 
to Pressure Ratio versus Time for an Infinite System; 
Pseudoplastic Fluids, d, = 2 
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derivative plots also equals v. In fact, the late-time value of the wellbore pressure-derivative 
function equals the product of the wellbore pressure function and the parameter v. And 
therefore, a plot of dln(p,,p)/din(tp) against tp would asymptotically tend to a constant value 
of v at late times (see Figure 10-4). An identical observation was also made in the study of 


Chang and Yortsos (1990) for the case of Newtonian fluid flow in a fractal object. 


There are certain advantages in using the pressure derivative group, d/n(py,p)/dln(tp), 
when constructing type curves for various flow situations (Duong, 1989). Duong has 
shown that by combining the pressure and pressure-derivative functions, a single set of type 
curves may be constructed by using the pressure/pressure-derivative ratio (PDR). By 
plotting dimensionless PDR (vertical axis) versus dimensionless time (horizontal axis), the 
vertical scales for both type-curve and field data are found to be identical. This alignment of 
the vertical scale constrains the type-curve match on the vertical axis. In this study, use is 
made of the dimensionless pressure-derivative/pressure ratio, which also serves the purpose 
of automatic alignment of one scale of the type curves and field data plots. This can be 
demonstrated by considering Equations (4-38) and (4-40) which can be rewritten, 


respectively, as 
Pwo = a(Pi- Py) = a(Ap) (10-1) 
tp = bt (10-2) 


where a and b are constants. The pressure-derivative, then, is given as 


tpDyp = atdp’, (10-3) 
and thus, 
din(Pwp) _ Apt (10-4) 


d \n(tp) Ap 
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Figure 10-5: Dimensionless Wellbore Pressure versus 
Time in an Infinite System; Dilatant Fluid, dy = 1.75 
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Figure 10-5 demonstrates the transient wellbore pressure behaviour for different 
values of n (dilatant fluid) and for d, = 1.75. Figure 10-5 shows that the dimensionless 
pressure decreases with increasing values of n, at any given time. Here also, the late-time 
linear behaviour in the log-log plot of pyp versus tp is apparent. Figures 10-6 and 10-7 
demonstrate the variation of pyp with time for various values d, of and for n = 0.75 and 
1.25, respectively. These figures show, as was shown by the previous figures, that with 
increasing magnitudes of v, the dimensionless pressure drop at a given time, increases. For 
example, in Figure 10-6, the lowermost curve corresponds to v = 0.111 and the uppermost 
curve to v = 0.256. Similarly, in Figure 10-7, the uppermost and the lowermost curves 


correspond to v = 0.2 and 0, respectively. 


In their groundbreaking study on mathematical modeling of single-phase transient 
flow of slightly compressible Newtonian fluids in fractal reservoirs, Chang and Yortsos 
(1990) have noted that with the sole exception of 2D cylndrical flow systems, the 
asymptotic pressure behaviour exhibited by any pressure-transient response model is of the 
power-law type DPwp ~ tp. The Warren and Root (1963) type double-porosity system 
exhibits the well-known asymptotic behaviour given by Dyp ~ /n tp. On the other hand, the 
single-fracture pressure response, both at early (linear flow period) and later (bilinear flow 
period) times, shows behaviours consistent with the description of pywp ~ th. The fractal 
reservoir model also exhibits such a linear behaviour in the late-time log-log plot of pressure 
versus time. It has been shown in the present study, as was demonstrated earlier by Ikoku 
and Ramey (1979), that the power-law description also applies for the case of non- 
Newtonian power-law fluid flow in a 2D cylindrical flow medium. For Newtonian flow in a 
fractal reservoir, the exponent v of the asymptotic pressure-time behaviour is only a function 
of the spectral dimension, d;. For non-Newtonian flow in a homogeneous reservoir, the 
exponent is only a function of the flow behaviour index, n. And has been shown earlier, for 


non-Newtonian flow in a fractal reservoir, the exponent v depends on both n and ds. 
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Figure 10-6: Dimensionless Wellbore Pressure versus 
Time for an Infinite System with varying dy; n = 0.75 
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Figure 10-7: Dimensionless Wellbore Pressure versus 
Time for the Flow of a Dilatant Fluid (n = 1.25) ina 
Fractal Reservoir; Infinite System, Varying dg 
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Figure 10-8 presents plots of pressure and pressure-derivative for two different cases 
of non-Newtonian fluid flow in a fractal medium. The upper pressure curve describes 
Newtonian flow in a fractal network and is characterized by a value of v of 0.125. The lower 
pressure curve describes non-Newtonian flow in a homogeneous medium and is defined by 
v = 0.111. The two lowermost curves in Figure 10-8 are the pressure-derivative plots and at 
late times, the parallel straight lines of pressure and pressure-derivative are separated by a 


distance equal to log(J/v). 


It can be inferred also from Figure 10-8, and from the other results presented so far in 
this chapter, that the transient pressure behaviour of non-Newtonian flow in a homogeneous 
medium and that of Newtonian flow in a fractal medium are similar in many ways. Thus, an 
important question that may arise in the context of the present study is how to separate the 
individual characteristics of the non-Newtonian fluid and of the fractal medium when the 
complexities of both are present in the response of the flow system. In other words, how is 
it possible to determine both n and ds from single-well test results if one can only calculate a 
value of the parameter v from the late-time pressure-time data? In order to find a possible 
solution to this question, one has to inspect the late-time behaviour of dimensionless rate for 


a constant-pressure wellbore condition. 


It may noted from Equations (5-30), (5-82) and (5-84) that the relationships between 


the constant rate solution (p,,,p) and the constant pressure solutions (¢p, Op), in Laplace 


space, are given as 


1 


(10-5) 
ap lend.) 


Pyp(l, n, a.) = 


and 
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—— Pressure: d, =1.75,n=1 
--- Pressure: d, = 2, n = 0.75 


PwD> dpwp/d In(tp) 


nnn « Pressure-Derivative: d, = 1.75,n=1 
-'«= Pressure-Derivative: d, = 2, n = 0.75 


Figure 10-8: Pressure and Pressure-Derivative 
Solutions for Two Different Cases 
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It may be interesting to note that, for an infinite outer boundary condition, these relationships 
would also hold for other flow situations considered in this study (see, for example, 
Equations (8-40), (8-53) and (8-54) for the case of a composite system; also see Equations 
(9-31), (9-44) and (9-46) for flow in the fracture/matrix system). Identical observations 
were made also in the classical study of van Everdingen and Hurst (1949) for the case of 


transient flow in a homogeneous flow medium. 


Figures 10-9 and 10-10 illustrate some transient dimensionless rate solutions. 
Equation (5-82) defines the dimensionless rate, gp, in Laplace space for the transient flow of 
a non-Newtonian power-law fluid in an infinite fractal medium. Equation (5-82) has been 
inverted to real space numerically by using the Stehfest algorithm (1970) in order to generate 
data for Figures 10-9 and 10-10. Figures 10-9 and 10-10 show that at early times the 
dimensionless rate curves are merged together, as predicted by Equation (5-85). With 
increasing time, the curves start diverging from each other with the degree of divergence 
depending on the values of n and d,. The value of dimensionless rate is higher at any given 
time for smaller magnitudes of v. For example, the uppermost curve in Figure 10-9 
corresponds to v = 0.111, and the lowermost one to v = 0.294. Similarly, the uppermost 
curve in Figure 10-10 corresponds to v = 0.125, and the lowermost curve to v = 0.263. At 
late times, the rate solutions exhibit linear behaviour in the log-log plot, as predicted by 
Equation (5-87). Thus, if the late-time rate versus time data is available, then it is possible to 
calculate the slope of this linear segment. From (5-87), the magnitude of the slope, for a plot 
of dimensionless rate (qp) versus dimensionless time (tp), is given by v. However, by 
combining Equations (5-80) and (5-87) it is clear that, a late-time log-log plot of 


dimensional rate (i.e., g) versus dimensional time (t) would be characterized by a straight 
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Figure 10-9: Dimensionless Rate versus Time in an 
Infinite System; n = 0.75 
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Figure 10-10: Dimensionless Rate versus Time 
in an Infinite System; dg = 1.75 
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line, the slope of which has a magnitude of v/n. Thus, it may be possible to calculate the 
values of both n and d,, provided both late-time rate (in a constant-pressure situation) and 


pressure (in a constant-rate situation) data are available. 


So far in the present chapter, various aspects of the late-time behaviour of transient 
pressure and rate solutions have been discussed. It has been demonstrated that, at relatively 
large times, the dimensionless wellbore pressure solution can be approximated by Equation 
(5-47). However, at fairly short times, the Laplace space wellbore pressure solution (given 
by Equation (5-30)) can be approximated very well by Equation (5-43), which is expressed 
in real space. In order to demonstrate the close approximation of Equation (5-30) by 
Equation (5-43), Figures 10-11 and 10-12, exhibiting comparisons of py values calculated 
from these two equations, are presented. It is clear that there is good agreement between the 
two solutions for tp greater than 100 and for v less than or equal to about 0.5. For larger 
values of v, the time when the two solutions match increases. It may be noted that a value of 
tp > 100 translates to real time of the order of a few minutes. Thus, for all practical 
purposes, the wellbore pressure solution can be represented by Equation (5-43), so long as 


the condition v < 0.5 is met. 


It is also of importance to be able to determine the times at the beginning of the log-log 
straight lines in the wellbore pressure versus time plots. Various types of correlations have 
been attempted in order to define the dimensionless time, at which the log-log straight lines 
appear, as a function of parameter v only. It has been found in this study that, in order to 
define a proper correlation of the said type, the dimensionless wellbore pressure and time 


should be redefined as follows: 


D = Ppyp(n+1—-nd,) (10-7) 


pettainti- (1d, (10-8) 
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—— Equation (5-30) 


- = - Equation (5-43) 


n=], G--1, v=0.5 


n=0.75, d,=2, v=0.111 


Figure 10-11: Comparison of Wellbore Pressure 
Solutions Using Equations (5-30) and (5-43) 
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=== Equation (5-30) 
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Figure 10-12: Comparison of Wellbore Pressure 
Solutions Using Equations (5-30) and (5-43) 
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Then the time (¢*), at which the pressure solution (p*) is within 2% of the log-log straight 
line behaviour, may be expressed as a function of parameter v, as shown in Figure 10-13. It 
should be noted that the correlation exhibited in Figure 10-13 is valid only for pseudoplastic 
fluids and for 1 < d, < 2. An explanation for why Equations (10-7) and (10-8) should be 
used in order to obtain the shown correlation may also be given by examining the large-time 
approximation of the wellbore pressure solution (Equation (5-43)). From Equation (5-43) it 


is clear that, after relatively short times, the wellbore pressure solution may be expressed as 


peas 10-9 
A ay ee eee 


when the wellbore pressure and time are expressed by Equations (10-7) and (10-8), 


respectively. 
10.1.2 Comparison of Analytical and Numerical Solutions 


In Chapter VI, a finite-difference scheme was presented for solution of the nonlinear 
form of the partial differential equation governing the transient flow of a non-Newtonian 
power-law fluid in a fractal reservoir. The main reason for attempting a numerical solution 
of the nonlinear partial differential equation, Equation (6-2), is to explore the consequences 
of the approximation given by Equation (4-29). This approximation resulted in the analytical 
solution for dimensionless wellbore pressure solution, given by Equation (5-30). The 
approximation introduced some incompressibility into the system (Ikoku, 1978; Ikoku and 
Ramey, 1982) and also removed the dependence of the analytical solution on the fractal 
dimension, dy, which may result in an inaccurate theoretical prediction of the model 
behaviour. Thus, the effects of this approximation should be analyzed by a direct 
comparison of the analytical and numerical solutions for transient dimensionless wellbore 


pressure during the flow of a power-law fluid in a fractal medium. 
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Figure 10-13: Plot of Dimensionless Time at the 
Beginning of Log-Log Straight Line Against 
Parameter v 
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Attempts were made to solve Equation (6-2) for different values of the flow behaviour 
index, n, and the spectral dimension, d,. The effect of the fractal dimension, d,, on the 
numerical solutions was also studied. Numerical solutions could not be obtained for values 
of n smaller than 0.3 and values of d, smaller than 1. In fact, the early-time numerical results 
for n values of 0.3 and 0.4 were not meaningful. The results presented in this section are 
based upon the following ranges of values of the parameters: 0.5 < n < 1, d, => 1 and 1.75 

< dy S 2. For the numerical solutions, the spatial increment was maintained at Ax = 0.1 


and a time increment of 10 was taken between printed results. 


Figure 10-14 presents a comparison of numerical and analytical solutions for n = 0.75, 
d= 2 and for three different values of d;. The curves (from top to bottom) are characterized 
by d, = 1, 1.5 and 2, which correspond to values of v = 0.5, 0.294 and 0.111, respectively. It 
is clear from these curves that, at large times, the difference between the analytical and 
numerical solutions is small. Moreover, the error in the analytical solution is seen to 
decrease, at any time, with an increasing magnitude of d,. Figure 10-15 demonstrates the 
temporal variation of the difference between the analytical and numerical solutions with 
varying n. The values of dy and d, are 2 and 1.75, respectively, for all the curves displayed in 
Figure 10-15. The graphs (from top to bottom) are characterized by n = 0.6, 0.8 and 1.25, 
which correspond to values of v of 0.239, 0.186 and 0.034, respectively. Here also, as in 
Figure 10-14, the error in the analytical solution decreases gradually with increasing time. 
The analytical wellbore pressure solution, given by Equation (5-30) for an infinite system, 
does not show any dependence on dy. In order to study the variation of pwn, evaluated 
numerically, with d,, Figure 10-16 is presented. The analytical solution for n = 0.5 and d; = 
2 is compared with three different numerical solutions defined also by dy values OFZ 71.85 
and 1.75. As Figure 10-16 shows, the three different numerical solutions can be 


distinguished in the log-log plot only upto a dimensionless time of about 50; beyond this 
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Figure 10-14: Comparison of Analytical and Numerical 
Solutions for Dimensionless Wellbore Pressure Variation 
(n=. 0075. d-=—2) 
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Figure 10-15: Comparison of Analytical and Numerical 
Solutions for Dimensionless Wellbore Pressure Variation 
(Ofona0c— 1:75) 
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Figure 10-16: Comparison of Analytical and Numerical 
Solutions for Dimensionless Wellbore Pressure Variation 
(n = 0.5, dg = 2) 
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time, the numerical solutions appear as one curve, and the analytical and numerical solutions 
are more or less parallel to each other. Figure 10-17 presents two different sets of analytical 
and numerical solutions; as expected, the upper curve corresponds to a higher value of v (= 
0.263) and the lower one to a smaller value of v = 0.006. Here also, the difference between 
the analytical and numerical solutions is larger at smaller times and for larger values of v. 
Figure 10-18 presents another graph describing the effect of n on the difference between the 
two solutions. Here, the values of n were chosen such that they are fairly close to that for a 
Newtonian fluid; also, dy= d, = 2. Figure 10-18 clearly demonstrates that the difference 
between the analytical and numerical solutions is practically negligible for tp > 100 for both 


sets of solutions. 


It may be observed, from Figures 10-14 through 10-18, that after large enough times 
(tp > 200), the numerical and analytical solutions appear to possess the same shape. Also, 
the difference between these two solutions generally seems to decrease with increasing times 
and decreasing magnitudes of the parameter v. The results from these graphs may be 
summarized by comparing the difference between the numerical and analytical solutions at 
tp values of 102 and 104 for all the curves presented in Figures 10-14 through 10-18. Such a 
comparison is presented in Table 10-1, where the relative difference is defined as the ratio of 
the magnitude of the difference between the numerical and analytical solutions to the value 
given by the numerical solution. As Table 10-1 shows, for a pseudoplastic fluid the relative 
difference is small at large times, provided the magnitude of n is large (i.e., close to 1, the 
Newtonian limit). For example, for n = 0.5 and df= d, = 2, the relative difference is about 
10% at tp = 104. Keeping everything else the same, the relative difference is seen to be 3% 
for n = 0.75 and < 1% for n = 0.9. Also, for n > 1, the relative difference is observed to be 
small at large times, provided the magnitude of v is small. Again, the analytical solution 
appears to yield good results if d, is large, every other parameter being the same. Finally, the 


large-time relative difference seems to be relatively insensitive to small variations in df , at 
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Figure 10-17: Comparison of Analytical and Numerical 
Solutions for Dimensionless Wellbore Pressure Variation 
(df = 2) 
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Figure 10-18: Comparison of Analytical and Numerical 
Solutions for Dimensionless Wellbore Pressure Variation 
(df= dg = 2) 


* 
i 
re i= 
} = >, 
“ >, 
i 
{ 
bs] 
ed Ee ‘. 
J 
*, 
: . be 
Sy 
x, 
», 
4 
ra as | 7 - dae A ae ee | : oan — ss wn -_ aos 


loko fag tested Tie tea BE-OL 
avitere fy Misa melee cel ego HT mab 
a te = b) 


_ 


Table 10-1: Comparison of Analytical and Numerical Solutions 


at Two Different Dimensionless Times 


Relative Difference, %, @ 


dr d, Vv tp = 102 tp = 104 
2.00 1.00 0.500 2 9 
2.00 1.50 0.294 6 =) 
2.00 2.00 Osh 4 3 
2.00 15 i239 10 8 
2.00 es: 0.186 4 3 
2.00 gis 0.034 1 1 
2.00 2.00 0.200 12 10 
1.85 2.00 0.200 12 10 
be7D 2.00 0.200 12 10 
2.00 175 0.263 14 10 
2.00 1.66 0.006 # Z 
2.00 2.00 0.048 1 << ih 
2.00 2.00 0.005 1 =205 
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least within the range of values of the various parameters studied. Thus, so long as the 
magnitude of n is larger than 0.75, that of dr is greater than 1.75 and that of d, is larger than 
1.50, the relative difference between the large-time analytical and numerical solutions would 
not be more than approximately 5%. More importantly, under these conditions, the relative 


difference between the slopes of the large-time analytical and numerical solutions would be 


even smaller. 
10.1.3 Analytical Solutions: Finite Fractal Medium 


In this section, some of the results obtained for the case of a finite fractal flow medium 
are presented. Two different situations are considered here: the first case involves a bounded 
circular fractal reservoir with both closed and constant pressure outer boundaries, and the 
second case involves a two-zone composite reservoir situation with the inner zone being a 
fractal medium and the outer zone a homogeneous region. The results presented in this 


section pertain to a constant flow rate condition imposed at the wellbore. 


Figure 10-19 is a graph of dimensionless pressure and pressure derivative for the case 
of a closed circular reservoir with n = 1, d, = 2 and rgp = 100. Three different plots of 
pressure and pressure derivative are displayed in Figure 10-19 for dp= 1.5, 1.75 and 2. The 
corresponding plots of d In(p,,p)/d In(tp) are displayed in Figure 10-20. Figure 10-19 
shows that at early times, in the infinite-acting stage, the pressure behaviour is independent 
of dy. With increasing time, the effect of the finite outer boundary is felt by the pressure 
transients; it is interesting to note that the smaller the value of dp, the earlier is the time when 
the pressure transients respond to the boundary effects. A similar behaviour is also exhibited 
by the pressure derivative, which attains the same value as the pressure response at large 
times when pseudosteady-state is attained (as shown by the unit-slope line). This is also 


shown in the late-time behaviour of d In(p,,p)/d In(tp). 
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—— Pressure Solution 


== Pressure-Derivative Solution 


Figure 10-19: Plots of Pressure Drop and Pressure- 
Derivative for a Closed Circular System: n = 1, 
ds = 2s reD = 100 
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Figure 10-20: Effect of dp on din(pwp)/dIn(tp) 
for a Closed Reservoir 


127 


Pied 7 
:: - iy 
ooo 7 — 
j > - ote ; | 
| ' z 
. vit i %S 
ae a. 
: i ve 
} ; Ta 
| > f : 
» 7. 
, a 


er 
~ 
7 


ak 


———— 
a ; 


Figure 10-21 presents a graph of dimensionless pressure and pressure derivative for a 
closed reservoir with n = 1, d, = 2 and r,p = 200. Here also, the plots of pressure and 
pressure derivative are characterized by dy values of 1.5, 1.75 and 2. The corresponding plots 
of dIn(p,p)/d In(tp) are shown in Figure 10-22. It may be seen from Figures 10-19 and 10- 
21 that the plots of pressure and pressure derivative for dg = 2 and rep = 100 and those for 
d¢ = 1.75 and rep = 200 are almost identical. In other words, for dy < 2, the pressure 


response of the system may be misinterpreted to represent a smaller reservoir size if one 


assumes dy to be equal to 2. 


Figure 10-23 presents plots of pressure, pressure derivative and d In(py,p)/d In(tp) 
solutions for a constant pressure outer boundary condition. The plots correspond to n = 1, d, 
= 2 and rep = 200 and dy = 1.5, 1.75 and 2. Here, at large times (tp > 105), the effect of dy 
on the pressure response is similar as in the case of a closed reservoir. For the same value of 
rep, the large-time pressure solution is smaller for a smaller value of dy. This may result in 
an underestimation of reservoir size if one assumes d- = 2, when in reality dris less than 2. 
The large-time difference between the three plots is displayed more clearly by the pressure 
derivative plots which, however, show the same general trend of a sharp decrease with 


increasing time due to the gradual attainment of steady-state. 


In the present study, the wellbore pressure solutions for the case of a composite 
reservoir will be displayed mainly by the pressure derivative response of the system. 
Because of its detail-enhancement capabilities, the pressure derivative has been used to a 
large extent in recent studies to analyze the behaviour of composite reservoirs of various 
geometries (Ambastha and Ramey, 1989; Stanislav et al., 1992). The results presented in 


this section are directed at studying the sensitivity of the response of the system to 
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—— Pressure Solution 


‘== Pressure-Derivative Solution 


Figure 10-21: Plots of Pressure Drop and Pressure— 
Derivative for a Closed Reservoir: n = 1, d, = 2, 
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Figure 10-22: Effect of dg on din(pwp)/dIn(tp) 
for a Closed Circular Reservoir 
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Figure 10-23: Effect of dpon the Pressure Behaviour 


of a Circular Reservoir with Constant Pressure 
Outer Boundary: n = 1, dg = 2, rep = 200 
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parameters such as the dimensionless size of the inner region (a), the permeability ratio 


(K,/k’), ratio of permeability-to-porosity-ratio (F = “oer dy, and ds. 

The variation of dimensionless wellbore pressure with time is shown in Figure 10-24 
for a two-zone composite system. The log-log plot is generated to illustrate the effects of 
permeability ratio and parameter F on the system response. Figure 10-24 shows that for a 
fixed value of permeability ratio, dimensionless pressure is higher at any time for higher 
values of the porosity ratio (@,,/6), even though the plots have similar shapes at relatively 
large times. Furthermore, for a fixed value of porosity ratio, the dimensionless pressure is 
significantly higher for higher values of the permeability ratio. The effects of permeability 
ratio and F on the system behaviour may be better illustrated by the derivative response, as 


shown in Figure 10-25. The following is apparent from Figure 10-25: 


i) After the inner region flow is terminated, transitional flow occurs when the derivative 
goes through a maximum corresponding to the values of permeability ratio and F; 
subsequently, the derivative declines and, at late times, it demonstrates the constant 


behaviour characteristic of an infinite homogeneous medium. 


ii) The permeability ratio has a strong impact on the time when the maximum in the 


derivative occurs and also on the magnitude of the maximum derivative. 
iii) The transition period seems to end sooner for the low permeability ratio case. 


iv) The parameter F has a mild influence on when the maximum in the derivative 


occurs, but has a stronger influence on the value of the maximum derivative. 
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Figure 10-24: Effect of Permeability and Diffusivity 
Ratios on Pressure Solution for a Composite 
Svstemians=1l, d6='2, de =e sad 
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Similar behaviour of derivative curves has also been observed for composite homogeneous 


reservoirs of cylindrical (Ambastha and Ramey, 1989) and elliptical (Stanislav et al., 1992) 


flow geometries. 


Figures 10-24 and 10-25 correspond to n = 1, de=2, d, = 1.5 and a = 2. Figure 10-26 
presents a log-log plot of pressure derivative versus time identical in every respect to Figure 
10-25, except that the former corresponds to a = 10. Comparing these two figures, it may be 
observed that even though the magnitude of the maximum derivative (during the transitional 
period) and the late-time constant value of the derivative are not affected by the value of a, 
the time when the maximum derivative is reached and that when the flattening out of the 
derivative takes place are dependent upon the magnitude of a. As expected, for larger values 
of a, the maximum in the derivative is reached later; also, for larger magnitudes of inner 


zone size, the time when the flattening out of the derivative begins is seen to increase. 


Another aspect of the derivative plots that is of interest is its early-time behaviour. It 
has been shown by Ambastha and Ramey (1989) that, in the absence of wellbore storage 
and for a sufficiently large inner-zone size, an infinite-acting behaviour is shown by the 
derivative at early times (a flat derivative curve); obviously, the magnitude of this constant 
derivative value does not depend upon the permeability and porosity ratios, size of the inner 
zone, and so forth. For the case of an elliptical system, the early-time derivative plot would 
not consist of a zero-slope line, because of the non-cylindrical nature of flow at early times 
(Stanislav et al., 1992). It was not possible in this study to determine conclusively the 
behaviour of the derivative curves at very early times (tp < 1), because of some numerical 
instabilities in the calculation of the pressure derivative from the Laplace space solution at 
small times. However, it is not difficult to observe from Figures 10-25 and 10-26 that, at 
early times, the pressure derivative solutions would not exhibit a zero-slope derivative line; 


this is not difficult to visualize based on the discussions presented in previous sections 
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Figure 10-26: Effect of Permeability and Diffusivity 
Ratios on Pressure Derivative for a Composite 
System: n = 1, dr= 2, ds = 1.5, a = 10 
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regarding the pressure derivative behaviour for a fractal medium. However, these early-time 
characteristics of pressure derivative are of little practical importance because of the presence 
of wellbore storage effects, which may obscure any early-time telltale "signature", as it were, 
of the system. On the other hand, at large enough times, the theoretical pressure derivative 
behaviour for different types of composite reservoirs are quite similar, as has been 
mentioned earlier. Thus, without a proper foreknowledge of the nature of the composite 
system, it is quite possible to match the field data with results from a model that does not 


represent the actual system at all. 


Figures 10-24 through 10-26 were generated by assuming constant values of dy and 
d,. Figures 10-27 and 10-28 present the effects of dy and d;, respectively, on the pressure 
derivative behaviour of a composite system. Figure 10-27 shows that a smaller value of d¢ 
may make the composite system appear as one with a smaller magnitude of permeability 
ratio, if it is assumed that d= 2. Similarly, it may be seen from Figure 10-28 that a smaller 
magnitude of d, may be misinterpreted as the system having a higher permeability ratio, if it 


is incorrectly assumed that d, = 2. 


It has been shown in the previous graphs that the pressure derivative of a composite 
system depends upon a number of parameters. It was proposed in previous studies on 
pressure transient analysis of composite reservoirs (Ambastha and Ramey 1989; Stanislav 
et al., 1992) that the time coordinate may be redefined as tap = tp/a’, in order to eliminate 
the dependence of the pressure behaviour on the size of the inner region. In this study also, 
the dimensionless time is redefined in a similar way, in order to investigate the possibility of 
removing the dependence of the pressure derivative on a. Figures 10-29 and 10-30 present 
derivative curves plotted against trp. In Figure 10-29, the values of dy and d, are the same, 
viz. 2. It is clear from Figure 10-29 that the correlation holds quite well for relatively large 


magnitudes of a. A similar observation was made also by Ambastha and Ramey (1989). 
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Figure 10-27: Effect of dp on Pressure Derivative 


for a Composite System: n = 1, ds = 1.5, 
a= 55 Ky/k' = 100 
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Figure 10-28: Effect of d, on Pressure Derivative 
for a Composite System: n = 1, dg=2,a=5, 
esd ae 100 
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Figure 10-29: Pressure Derivative Behaviour for a 
Composite Reservoir: Different Inner-Zone Sizes 
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Figure 10-30 shows that for dy= 1.75 and d, = 1.25, the correlation holds for a = 2 and 5 
and not for a = 10. Again, for dp = d, = 1.5, the correlation does not seem to work for the 
three values of a, at least for trp < 104. Thus, it is reasonable to conclude from Figures 10- 
29 and 10-30 that a correlation of pressure derivative versus tgp does not work for all 
combinations of dr and d,, and alternative correlations must be sought for such composite 


reservoirs. 
10.2 Transient Flow Through a Fractal Fracture/Homogeneous Matrix System 


The previous sections dealt with the case of the matrix not participating in the system 
pressure response. This section presents a discussion of results obtained for the case of a 
fracture/matrix system (i.e., the matrix has connected porosity and is in pressure 
communication with the fracture network) with the fracture network being characterized by 
dg =2 and d, < 2 and the matrix being homogeneous. The results of calculation of the 
dimensionless wellbore pressure for a constant rate condition at the wellbore and an infinite 
reservoir could be used to analyze the characteristics of pressure drawdown curves in such a 
double-porosity reservoir for single-well situations. Figure 10-31 presents a basic pattern of 
such pressure drawdown curves for an infinite system on a log-log plot. As Figure 10-31 
shows, a complete curve may generally be represented by three different flow regimes or 
intervals. The early-time flow regime (interval I), characterized by the log-log straight line, 
exhibits a rapid drawdown in the fractal fracture system. Obviously, for flow only in the 
fracture network, the drawdown curve would maintain the same trend of interval I. 
However, for a fracture/matrix system, the transitional stage (shown by interval II) appears 
at the end of the first flow regime. The transitional stage involves a slower pressure 
drawdown in the fracture system because of increasing participation of the tighter matrix in 
the flow process. Interval III, the late-time straight line, develops after an equilibrium is 


reached between the fracture and matrix pressures. The straight line of interval I is parallel 
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Figure 10-31: Pressure Drawdown Curves 
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to that of interval I because of the assumption of df =i Nnleye df #2, the slopes of these two 
straight lines would be different, even though the difference would be small (Chang and 
Yortsos, 1990). For a sufficiently large distance to the outer boundary, the pressure response 


would show the boundary effects after interval III. 


Figure 10-32 presents sensitivity runs carried out to determine the effects of parameter 
v on the transient dimensionless pressure response. For the set of pressure curves presented 
in Figure 10-32, the values of @ and 4 are set to be 0.01 and 10 -5, respectively. The 
characteristics of the "double-porosity" system are exhibited more clearly by the pressure 
derivative curves, as shown in Figure 10-33. The derivative curves (for v 0) show the 
typical power-law dependence on time of the late-time segments, as is expected for fractal 
systems. The well known "V" shaped derivative behaviour for the Warren and Root (1963) 
type double porosity system (v = 0) appears to be quite similar to that for the fractal 
systems. However, the late-time pressure derivative for a homogeneous system stabilizes 
asymptotically to a finite value, unlike that for a fractal system. Figure 10-34 presents plots 
of the pressure derivative group, d In(py,p)/d In(t), with values of v ranging from 0 to 0.5. 
Clearly, for large magnitudes of v, the values of the derivative group would stabilize at 
constant values of v, unlike the case of a homogeneous fracture/matrix system, where the 
derivative group would decrease indefinitely with increasing times. Identical results have 


been presented previously, by numerical means, in the study of Chang and Yortsos (1990). 


Figure 10-35 presents the sensitivity of the wellbore pressure solutions to @, with @ 
ranging from 0.001 (relatively small amount of storage in the fracture system) to 1 (flow 
only in the fractures). Figure 10-35 shows that the early- and late-time pressure curves have 
the same slopes irrespective of the magnitude of @, the parameter w, however, has a strong 


influence on the duration of the transition period. As may be expected, the smaller the 
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Figure 10-32: Plot of Pressure Drop for a Fracture 
/Matrix System: @ = 0.01, A = ter 
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Figure 10-33: Plot of Pressure-Derivative for a 
Fracture/Matrix System: w = 0.01, A = ners 
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Figure 10-34: Plot of d In(pwp)/d In(t) for a 
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magnitude of @ relative to one, the earlier is the time when the transition stage sets in (so 
that the longer is the duration of the transition period) and the larger is the pressure drop 
between the early- and late-time segments. The effect of w on the pressure derivative is 
shown in Figure 10-36, for v = 0.25 and A = 10 -5. The transitional stage (when the pressure 
response shows a period of gradual decrease and subsequent increase in slope) appears as a 
dip in the derivative curve; the smaller the magnitude of @, the sharper is the minimum and 


the smaller is the value of t when the dip in the derivative takes place. 


Figures 10-35 and 10-36 demonstrated the variation of the pressure and derivative 
responses, respectively, against T, where thas been defined as tT = @fp. Figures 10-37 and 
10-38 are presented to study the effects of @, when the pressure and derivative responses are 
studied against dimensionless time tp. Figure 10-37 shows that at early times, the curves for 
different values of w merge together, as, during this period, flow takes place only through 
the fracture system. At late times, the pressure curve for the fracture/matrix system in the 
log-log plot is parallel to the hypothetical straight line for @ = 1 (shown by the dotted line). 
This is strictly true only for dr= 2, and for such a situation, the vertical distance between the 
two asymptotes depends only on two parameters, namely, v and @. However, when d; #2, 
the magnitude of the pressure drop not only depends on other parameters but is also time- 
dependent; therefore, when dr #2, it may not be possible to obtain a reliable estimate of the 
magnitude of @ from an observed value of the vertical pressure drop. On the other hand, 
when dy = 2, the late-time equations for wellbore pressure (Equation (5-47) for flow only in 
the fracture network and Equation (9-37) for flow in the fracture/matrix system) may be 
analyzed to show that the magnitude of the pressure drop varies with vlog(1/@). Thus, 
knowing the magnitude of v and that of the vertical pressure drop, the value of @ may be 


calculated. 
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Figure 10-36: Effect of @ on dpwp/dIn(t) for Fracture 


/Matrix System: v = 0.25,2 = 10> 
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Figure 10-37: Effect of @ on Plot of pywp vs. tp 


for Fracture/Matrix System: v = 0.25, \ = ion 
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Figure 10-38 shows the effect of won the pressure derivative behaviour. It may be 
noted from Figure 10-38 that the parameter @ has a weak influence on the dimensionless 
time when the pressure derivative reaches its minimum. This time, however, is strongly 
influenced by the parameter A, as is shown in Figures 10-39 and 10-40. The parameter A 
indicates the strength of the exchange rate between the matrix and the fracture system 
(Chang and Yortsos, 1990) and, therefore, the smaller the magnitude of A, the later is the 
onset of the transitional stage. It has been found in the present study, that within the ranges 
of the various parameters studied, the time when the wellbore pressure derivative reaches its 
minimum during flow in a fracture/matrix system depends on A and, to some extent, on @. 


The following empirical relationship is found to hold between this time and the parameters A 


and @: tp = In(1/@)/A. 
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Figure 10-39: Effect of A on Pressure Solution 
for a Fracture/Matrix System: w = 0.01, v = 0.25 
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Figure 10-40: Effect of on Pressure-Derivative 
for a Fracture/Matrix System: v = 0.25, m = 0.01 
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Chapter XI 


CONCLUSIONS 


A mathematical model describing the single-phase transient flow of a slightly 
compressible, non-Newtonian, power-law fluid in a fractal fracture network has been 
presented. The nonlinear partial differential equation governing such flows has been shown 
to reduce to the standard diffusivity equation as a limiting case. Approximate analytical 
solutions of the governing equation, with constant rate and constant pressure inner boundary 
conditions, have been presented for infinite and finite systems and for a composite reservoir 
case. Analytical solutions have also been presented for a special case of the flow situation 
where the matrix participates in the flow process along with the fracture network. Numerical 
solutions of the rigorous flow equation were compared with the approximate analytical 


solutions for flow in an infinite fractal reservoir. 


Based on the results obtained in this study, the following conclusions can be drawn: 

1) log-log plots of wellbore pressure and pressure derivative may be used as diagnostic 
plots to identify either the fractal nature of a reservoir or the power-law characteristics of 
the flowing fluid; 

ii) for power-law flow in a fractal reservoir, the slopes of late-time straight lines in the log- 
log pressure-drawdown plots will depend on both the flow behaviour index of the power- 
law fluid (n) and the spectral dimension of the fractal reservoir (d,). It is necessary to 
analyze the late-time behaviour of both transient pressure and rate data to calculate the 
values of n and d, under such flow situations; 

iii) the approximate analytical solution for pressure-drawdown behaviour in an infinite 


fractal system is seen to be within acceptable error, when the parameters have the 


following ranges: 0.75 $< n<J1,175 S$ dy Seo .a0d Lops tare sae 
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iv) the large-time pressure response of a finite, areally heterogeneous system with a fractal 
structure may be misinterpreted to represent a smaller reservoir size if analyzed using 
methods derived for homogeneous systems; 

v) methods presented in previous studies for the analysis of composite reservoirs may 
result in gross misinterpretation of the nature of the flow system if the inner region of the 
composite system consists of a fractal medium. The analysis of drawdown data from such 
reservoirs with techniques that do not account for the fractal characteristics of the inner 
zone can lead to incorrect reservoir property estimates; 

vi) for Newtonian flow in a reservoir comprising a fractal fracture network and a 
homogeneous matrix, the spectral dimension, d,, and the interporosity flow parameters, @ 


and A, may be estimated provided dp= 2 and the magnitude of ois small. 


One can make many recommendations for further research in this area. The more 
important ones are: 

i) to validate the analytical models presented in this study for finite and composite 
reservoirs by comparing the analytical solutions with rigorous numerical results and to 
subsequently define the ranges of magnitudes of the various parameters within which the 
analytical solutions compare favourably with the numerical ones; 
ii) to investigate the pressure transient behaviour of a fractal fracture/homogeneous matrix 
system for non-Newtonian fluids using numerical techniques; 
iii) to develop extensions of the methods presented in this study to multiple well situations 
and interference testing. It has been shown (Chang and Yortsos, 1992) that approaches 
based on Green's functions may be appropriate for such purposes; 
iv) to develop an understanding of the relationship, if any, between theories developed for 


a fractional dimension medium (e.g., Barker, 1988; Doe, 1991) and those for a fractal 


object. 
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Appendix A 


COMPUTER PROGRAM FOR GENERATING DIMENSIONLESS PRESSURES 
FOR POWER-LAW FLOW IN AN INFINITE FRACTAL RESERVOIR 


*% THE PROGRAM LISTED HERE IS DEVELOPED x 
*% TO CALCULATE DIMENSIONLESS WELLBORE * 
x PRESSURE RESPONSE IN LAPLACE SPACE DURING x 
x THE PRODUCTION OF A NON-NEWTONIAN, POWER- x 
* -LAW TYPE FLUID AT A CONSTANT RATE FROM * x 
xx AN INFINITE RESERVOIR EXHIBITING FRACTAL * x 
** DISTRIBUTIONS OF POROSITY AND PERMEABILITY. ** 
xx K* 
* THIS PROGRAM COMPUTES THE DIMENSIONLESS ** 
* PRESSURE RESPONSE AND ITS DERIVATIVE USING ** 
* NUMERICAL LAPLACE TRANSFORM INVERSION WITH ** 
x THE STEHFEST ALGORITHM (STEHFEST, 1970). xx 
xx* xx 
KKK KK KK MAIN PROGRAM KKK KKK KK KK 


xxx kK DEFINITION OF VARIABLES *****x* xx 
KKK KK KKK KKK KKK KK KKK KKK KKK KKK KKK KKK KKEKKKKK KKK 


** INPUT VARIABLES ***** RR KKK KK KK KKK 


AN --> DIMENSIONLESS FLOW BEHAVIOUR INDEX (POWER- 
LAW PARAMETER) 

DS --> DIMENSIONLESS SPECTRAL DIMENSION 

Se-->)  LAPLACE) VARTABLE 

T(J) --> ARRAY OF DIMENSIONLESS TIME 

TD --> DIMENSIONLESS TIME 


xx OUTPUT VARIABLES KAEKKKKEKKKKKKKKKKK KKK KKK K 


PWD(J) --> VECTOR OF THE WELLBORE PRESSURE DROP 
DPWDL(J) --> VECTOR OF THE SLOPE OF THE WELLBORE 
PRESSURE RESPONSE 


KKK KKK KKK KK KK KK Kk ke kk kk kK kK Ke Kk KKK OK Kk KK KK 


IMPLICIT REAL*8(A-H, O-Z) 


DIMENSION TDST(10), TD(100), 
# PWD(100), DPWDL(100) 


COMMON /THE/AN, DS 


DATA TDST/1.00D0,1.50D0,2.00D0,2.50D0,3.00d0,4.000,5.0D0, 
# 6.00d0,7.0000,8.00d0/ 

M=777 

N=8 

READ (5,%*) AN,DS 

CALCULATE THE VALUES OF TD 

ITD=0 

DO 20 J=1,8 

DO 20 K=1,10 
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ITD=ITD+1 

IDI=TDST (K) * (10.0D0** (J-1) ) 

TD (ITD) =TDI 

CONTINUE 

NTD=ITD-1 

WRITE (6,30) AN,DS 

FORMAT ('The values of n =',1E9.4,2X,'and ds =',1E9.4) 
WRITE (6,40) 

FORMAT (3X, *tD", 5X, ‘PwbD", 5x, "dPWD/dintD") 


DO 50 J=1,NTD 

t=ANB) (3p) 

CALL INVERT (T,M,N,PD,PDP) 

PWD (J) =PD 

DPWDL (J) =T*PDP 

PRINT OUT COLUMNS OF TD, PWD AND DPWDL VALUES 


WRITE (6,60) T,PWD(J) ,DPWDL(J) 
BORMAT® (1E12°6;5X,1E12.6,5X,1E12.6) 
CONTINUE 
STOP 
END 

THE STEHFEST ALGORITHM 

KK KK KK RK KK KK RK KK RK IK Kk 

SUBROUTINE INVERT (TD,M,N,PD, PDP) 
THIS SUBROUTINE COMPUTES NUMERICALLY 
THE LAPLACE TRANSFROM INVERSE OF ANY FUNCTION 
TEN GAPLACE SPACE, E(S):- 
COURTESY OF: ANIL K. AMBASTHA (PET E 668; 1991) 


IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION G(50), H(50), V(25) 


NOW IF THE ARRAY V(I) WAS COMPUTED BEFORE, THE PROGRAM 
GOES DIRECTLY TO THE END OF THE SUBROUTINE TO CALCULATE 
F(S) - 


IF (N.EQ.M) GO TO 17 

M=N 
DLOGTW=0 . 6931471805599453D00 
NH=N/2 


THE FACTORIALS OF 1 TO N ARE CALCULATED INTO ARRAY G. 


G(1)=1.D0 

DO 1 I=2,N 
G(ty=G(I-1) *I 
CONTINUE 


TERMS WITH K ONLY ARE CALCULATED INTO ARRAY H. 


H(1)=2.D0/G(NH-1) 

DO 6 I=2,NH 

FI=I 

IF (I-NH) 4,5,6 

H (I) =FI**NH*G (2*1) / (G(NH-I) *G(I) *G(I-1) ) 
GO TO 6 

H (I) =FI**NH*G (2*I) / (G(I) *G(I-1)) 
CONTINUE 


THE TERMS (-1) **NH+1 ARE CALCULATED. 
FIRST THE TERM FOR I=1 
SN=2* (NH-NH/2*2) -1 
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THE REST OF THE SN'S ARE CALCULATED IN THE MAIN ROUTINE. 


THE ARRAY V(I) IS CALCULATED. 
DO 7 I=1,N 


FIRST SET V(I) =0 
V(I)=0.DO0 


THE LIMITS FOR K ARE ESTABLISHED. 
THE LOWER LIMIT IS K1=INTEG((I+1/2) ) 


K1=(I+1)/2 

THE UPPER LIMIT IS K2=MIN(I,N/2) 
K2=I 

TE e(KZ-NH) 6, 6; 9 

K2=NH 


THE SUMMATION TERM IN V(I) IS CALCULATED. 


DO 15 K=Ki,K2 

phate eK-T)e 12,13, 12 

ie Yea ey eles: ath 

V (I) =V (I) +H (K) / (G(I-K) *G(2*K-I) ) 
GOLTo. 15 

V (I) =V(1I) +H (K) /G (I-K) 

GO TO 15 

V (I) =V (I) +H (K) /G(2*K-I) 

CONTINUE 


THE V(I) ARRAY IS FINALLY CALCULATED BY WEIGHTING 
ACCORDING TO SN. 
V (I) =SN*V (T) 


THE TERM SN CHANGES ITS SIGN EACH INTERATION. 
SN=-SN 
CONTINUE 


THE NUMERICAL APPROXIMATION IS CALCULATED. 
PD=0 .D00 

PDP=0 .D00 

A=DLOGTW/TD 

DO 19 I=1,N 

ARG=A* I 

CALL LAP (ARG, PWDL, PDPL) 
PD=PD+V (1) *PWDL 
PDP=PDP+V (I) *PDPL 
CONTINUE 

PD=PD*A 

PDP=PDP*A 

RETURN 

END 


SUBROUTINE LAP HAS THE FUNCTIONS THAT ARE OF INTEREST; 
THE DIMENSIONLESS PRESSURE AND DIMENSIONLESS PRESSURE 
DERIVATIVE FUNCTIONS ARE EXPRESSED IN TERMS OF THE 
LAPLACE VARIABLE, S, WHICH CORRESPONDS TO VARIABLE 
"ARG" IN SUBROUTINE "INVERT". 


SUBROUTINE LAP (S,PWDL, PDPL) 

IMPLICIT REAL*8 (A-H,0-Z) 

COMMON /THE/AN, DS 
A=2.0D0*DSQRT (S) / (AN+1 .0D0-DS* (AN-1.0D0) ) 
V1=1.00d00-DS/ (AN+1.0D0-DS* (AN-1.0D0) ) 
V2=1.0D0-V1 
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67 


IF (V1 .EQ. 0.0D00) GoTo 66 


NN=1 

CALL DBSKES (V1,A,NN, BK1) 

CALL DBSKES (V2,A,NN, BK2) 
ANUM=BK1/DEXP (A) 
ADEN=BK2*S*DSORT (S) /DEXP (A) 

GOTO 67 

ANUM=DBSKOE (A) /DEXP (A) 
ADEN=DBSK1E (A) *S*DSQRT(S) /DEXP (A) 
PWDL=ANUM/ADEN 

PDPL=S*PWDL 


RETURN 
END 


167 


Ce OE @)1G)@) ug GG OlOrG) Ol@.@1@)@)@) | @) GO ai@ OO) G@:O) GO Oro @l@) O1@ OM) OLOlo ka Gi oleke 


Appendix B 


COMPUTER PROGRAM FOR NUMERICAL SOLUTION OF THE NONLINEAR 
PARTIAL DIFFERENTIAL EQUATION FOR FLOW OF A POWER-LAW 
FLUID IN AN INFINITE FRACTAL RESERVOIR 


** THE PROGRAM LISTED HERE USES THE DOUGLAS- *x 
x JONES PREDICTOR-CORRECTOR METHOD FOR THE x 
x NUMERICAL SOLUTION OF THE NONLINEAR * x 
x PARTIAL DIFFERENTIAL EQUATION GOVERNING ala 
* THE TRANSIENT FLOW OF A NON-NEWTONIAN wx 
«x POWER-LAW FLUID IN AN INFINITE FRACTAL ** 
** RESERVOIR. wx 
x* K* 
* AT TIME TD=0, THE DIMENSIONLESS PRESSURE * x 
x IS ZERO EVERYWHERE. THE FLOW RATE INTO THE ** 
ek FINITE-RADIUS WELLBORE IS CONSTANT, AND * 
** PRESSURE AT THE OUTER BOUNDARY IS ZERO AS * 
x THE DISTANCE TO THE OUTER BOUNDARY TENDS x 
x TO A VERY LARGE VALUE. he 
*x* x* 
x THE WELLBORE PRESSURE IS PRINTED AT EVERY * 
x TIME-STEP AFTER CONVERGENCE IS ACHIEVED AT ** 
** THE CORRECTOR STAGE. ** 
KKEKKKKKK MAIN PROGRAM RKKKKKKKKK KK 


xk KKKKKK DEFINITION OF VARIABLES ******* kx 
KKK KKK KKK KKK KKK KKK KKK KKK KE KKK KKKKKKKKK KKK K KKK 


** INPUT VARIABLES * **® RK KK RK KR RK I 


AN --> DIMENSIONLESS FLOW BEHAVIOUR INDEX (POWER- 
LAW PARAMETER) 

DELID) ——> SIZE OF THE LIME=STEP 

DF --> DIMENSIONLESS FRACTAL DIMENSION 

DS --> DIMENSIONLESS SPECTRAL DIMENSION 

DX --> SPATIAL INCREMENT 

EPSILON --> MAXIMUM ALLOWABLE RELATIVE DIFFERENCE IN PRESSURES 

FROM TWO DIFFERENT ITERATION LEVELS 

N --> NUMBER OF SPATIAL INTERVALS INTO WHICH THE SYSTEM 
IS DIVIDED BY THE GRID POINTS 

TD --> DIMENSIONLESS TIME 

TDMAX --> MAXIMUM VALUE OF DIMENSIONLESS TIME CONSIDERED 


**k OUTPUT VARIABLES (*% % 45% 4% 511% 1 9 I IK 
IC --> ITERATION COUNTER 


P(1) --> REQUIRED WELLBORE PRESSURE SOLUTION AT THE 


CURRENT VALUE OF TIME TD 
KK IK KK KK FO RO OK OK OK I I OR KIO IO II III kk 


IMPLICIT REAL*8(A-H, O-2Z) 


DIMENSION A(1001), B(1001), C(1001), D(1001), P(1001), 
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READ (5,*) N, DX, DELTD, TDMAX 
READ (5,*) AN, DF, DS, EPSILON 
WRIEE(67."(2 (A,G12..6))") 


WRITE) (6, "(2(A,G12.6)).*) 

me LDS 3=e" 7, Ds; ", EPSILON =)", EPSTLON 
eS Day 

Nal Sap al 

A3 = (DF/DS) * (AN + 1.0D0) / AN 

A2 = (DF/DS) ** ((AN + 1.0D0) / AN) 

AD = AN * (DF °= A3) 


DX2 = DX * DX 
ALPHA = 2.0D0 * DX2 / DELTD 


Cl = A3 * DX 

C2e= (ANU-71 20D0). /VAN 
C2Na=—=—C2 

C3"= (DF/DS) ** C2 

C4 = DxX* DF / DS 

C5 = A2 * ALPHA 

C67= Al ** Dx 


SETTING INITIAL AND BOUNDARY CONDITIONS... 
DODO et = 1, NL 

P(I) = 0.0D0 
PHALF (I) = 0.0D0 


PERFORM CALCULATIONS OVER SUCCESSIVE TIME-STEPS 
TD 0.0d00 
AUB, AWD) Ge" DANE AMD) 


NOW BEGINS THE PREDICTOR; THE COEFFICIENT ARRAYS A, B, C AND 
RIGHT-HAND SIDE VECTOR D ARE SET AND PRESSURES AT THE END OF 
HALF OF THE TIME-STEP, PHALF, ARE COMPUTED. 


B(1) = -(2.0D0+C3*C5) 
Di pies Ose C60) C3) *°C5 oP (1) —9 2 0D0* ca 
AiG) = —* 0). ODO 
Ctl) ZA0D0 
DO 60 I = 2, N 
Teo et 7 
TM1 = TF —<1 
AIM1 = IM1 
OP1 = (P(IM1) = P(IPL)) / (2.0D0*Dx) 


TPe(DPL” LE. 0.0D00) GO TO 30 

PD1 = 1.0D00 / DP1l 

B(I) = -(2.0D0+C5*DEXP (C1*AIM1) * (PD1**C2N) ) 

Dy eal DX2 SDP ie —4C54,* DERE (C1 *AIM)) =. (BDI *C2N) exer UL) 
GO TO 50 


CONTINUE 
IF (TD .EQ. DELTD) THEN 
B(I)=-2.0D0 


D(I)=0.0D0 
GO TO 50 
ENDIF 
NNo=2r = 1 


CAL] TRISOL(1, NN, A, B, C, D, PNEW) 
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DO 40 J = NN +1, N 
ENBW (0) © =°F (J) 

GO TO 70 

CONTINUE 


A(I) = 1.0D0 
Cel) = 1 0D0 
CONTINUE 
GUN aes 07 ODO 
CAPEIRISOL{ Nj2A,?B, 'C, D, 


DOPCOR TL = Fi N 
PHALF (I) = PNEW(1I) 


PNEW) 


THE CORRECTOR BEGINS NOW; IN THIS STAGE ALSO THE COEFFICIENT 
ARRAYS A, B AND C AND THE RIGHT-HAND SIDE VECTOR D ARE SET; 
SUBSEQUENTLY, THE PRESSURES AT THE END OF THE TIME-STEP 


ARE CALCULATED. 


Crag 
DL 
TODO EeeC4s * C6 — 4.0D0 * C4 
IEG, Se Ke Gp al 


DOelLs0 si -= 2, N 


f 
IN2IL GS ak Ge al 
Gk Sat al 
AIM1 = IM1 


DP2 = (PHALF(IM1) - PHALF (IP1) ) 


-2.0D0 * PHALF (2) +°2.0D0 * PHALF (1) 


PeetDr2. LE. 0.0D00) GO TO 100 


PD2 = 1.0D00 / DP2 


/ (2.0D0*DxX) 


B(Z) = -—(2.0D0+CS*DEXP (C1*AIM1) * (PDZ2**C2N) ) 
D(I) = -—(PHALF(IM1) + PHALF(IP1) - 2.0D00*PHALF (T) ) 


gel AIM)) © * AED2**C2N) * PT) 
GO TO 120 


CONTINUE 
NN eee 
CA11 TRISOL(1, NN, A, B, C, 


DO 110 J = NN +1, N 
PNEW(J) = PHALF (J) 


GO TO 140 


CONTINUE 
CONTINUE 


CAll) TRISOL (1, N, A,B, C, -D;, 


CHECK FOR CONVERGENCE 


ehCobeeCo,* P(h)  +e2. 


= (GS) ts DEAS ( 


+2. 0DUR* OB DAeae OrZ 


D, 


PNEW) 


PNEW) 


DIFF1 = (PNEW(1) - PHALF(1)) / PNEW(1) 
IF (DABS(DIFF1) .LE. EPSILON) GO TO 160 


Dor tou. 
PHALF (1) 


PNEW (TI) 


GO TO 90 
CONTINUE 


Thr 38 


BOSE 7 0.1 
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He) (PV 


P(I) = PNEW(I) 


PRINT CURRENT VALUE OF TIME, THE NUMBER OF ITERATIONS 
REQUIRED AND THE WELLBORE PRESSURE VALUE. 


WROTE OAs Gl2. G7 A la.) ) UTD (oun) DT eset Ce me 
WRITE 06,5 §(A,G1I276) Se P(E) 7 2) 

IP (TD. LE. TDMAX) GO TO 20 

STOP 

END 


SUBROUTINE TRISOL(I, N, A, B, C, D, PRESS) 

SUBROUTINE FOR SOLVING A SYSTEM OF LINEAR SIMULTANEOUS 

EQUATIONS WITH A TRIDIAGONAL COEFFICIENT MATRIX. 

IMPLICIT REAL*8 (A —- H,O = 2) 

DIMENSION A(*), B(*), C(*), D(*), PRESS(*), BETA(1001), 
GAMMA (1001) 

BEIA( I) =" B(I) 

GAMMA(I) = D(I) / BETA(TI) 

dhe Fe Ab Ge Tl 


DO" 10 K = IP1, N 
BETA (KJe= (B(K) —— A(R *-C{ Re) /* BETA(Ke—= A) 
GAMMA (K) = (D(K) - A(K) * GAMMA(K - 1)) / BETA(K) 
CONTINUE 
PRESS (N) = GAMMA (N) 
Mee Nae 
DO 20035 = 1, M 
Kee Neo) 


PRESS (K) = GAMMA(K) - C(K) * PRESS(K + 1) / BETA(K) 
CONTINUE 


RETURN 
END 
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Appendix C 


COMPUTER PROGRAM FOR GENERATING DIMENSIONLESS PRESSURES 
FOR POWER-LAW FLOW IN A FINITE FRACTAL RESERVOIR 


Pe THE PROGRAM LISTED HERE IS DEVELOPED ai 
ret TO CALCULATE DIMENSIONLESS WELLBORE ia 
“hg PRESSURE AND DERIVATIVE RESPONSE IN “ee 
a LAPLACE SPACE DURING THE PRODUCTION rales 
oaks OF A NON-NEWTONIAN, POWER--LAW TYPE a3 
ane FLUID AT A CONSTANT RATE FROM A FINITE rita 
3 FRACTAL RESERVOIR. Le! 
** x 
ests THIS PROGRAM COMPUTES THE DIMENSIONLESS ret 
rs PRESSURE RESPONSE AND ITS DERIVATIVE USING ‘stig’ 
aces NUMERICAL LAPLACE TRANSFORM INVERSION WITH rds 
Eats THE STEHFEST ALGORITHM (STEHFEST, 1970). tees 
*x* *x 
KKK KKK KK MAIN PROGRAM BOR OK KK kk 


kkk kK DEFINITION OF VARIABLES *********x 
KKK KKK KK KKK KKK KKK KKK KKK KKK KK KKKKKKEKKKKKKKKKK 


*k INPUT VARIABLES KRKKKKKKKKK KKK KKK KKK KKK KKK 


AN --> DIMENSIONLESS FLOW BEHAVIOUR INDEX (POWER- 
LAW PARAMETER) 
DF --> DIMENSIONLESS FRACTAL DIMENSION 
DS --> DIMENSIONLESS SPECTRAL DIMENSION 
RED --> DIMENSIONLESS RADIUS TO OUTER BOUNDARY 
S --> VARIABLE OF LAPLACE TRANSFORM W. R. T. TIME, TD 
T(J) --> ARRAY OF DIMENSIONLESS TIME 
TD --> DIMENSIONLESS TIME 


**k OUTPUT VARIABLES % ®% % 1K KKK IK I I II I IK 
PWD(J) --> VECTOR OF THE WELLBORE PRESSURE DROP 


DPWDL(J) --> VECTOR OF THE SLOPE OF THE WELLBORE 
PRESSURE RESPONSE 


KK KKK KK KKK KKK KKK KKK KKK KKK KKK KK KK KKK KKK KK KKK KKK KK KK KK 


IMPLICIT REAL*8 (A-H, O-Z) 


DIMENSION TDST(10), TD(200), 
# PWD(200), DPWDL(200) 


COMMON/THE/AN, DF, DS, RED 


DATA TDST/1.0D0,1.5D0,2.0D0,2.5D0,3.0D0,4.0D0,5.0D0,6.0D0, 
te ODO; 8. OD0/ 

M=777 

N=8 

READ (5,%*) AN,DF,DS,RED 
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60 
50 


CALCULATE THE VALUES OF TD 
ITD=0 

DO 10 J=3,10 

DO 10 K=1,10 

ITD=ITD+1 
TDI=TDST (K) * (10 .0D0** (J-1) ) 
TD (ITD) =TDI 


CONTINUE 

NTD=ITD-1 

WRITE (6,20) AN,DF 

FORMAT ('The value of n =',159.4,2X,'and df =',189.4) 


WRITE (6,30) DS,RED 

FORMATS (@ithe value of dse=!71Eb9. 4, 2x) land reD =", 18924) 
WRITE (6, 40) 

FORMAT "(SX,'tD', 3X, ‘PwD', 3X, *dPWD/dlntD") 


DO 50 J=1,NTD 

T=TD (J) 

CALL INVERT (T,M,N,PD,PDP) 
PWD (J) =PD 

DPWDL (J) =T*PDP 


PRINT OUT COLUMNS OF TD, PWD AND DPWDL VALUES 


WRITE (6,60) TD(J),PWD(J),DPWDL(J) 
FORMAT (1E12.6,3X,1E12.6,3X,1E12.6) 
CONTINUE 

STOP 

END 
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*** CLOSED OUTER BOUNDARY CASE .- RKKKKKKKKKKKKKKKK KK 


SUBROUTINE LAP HAS THE FUNCTIONS THAT ARE USED TO 
CALCULATE THE DIMENSIONLESS PRESSURE AND PRESSURE 
DERIVATIVE FUNCTIONS IN TERMS OF THE LAPLACE VARIABLE, 
S, AND THE SYSTEM PARAMETERS. 


SUBROUTINE LAP (S,PWDL, PDPL) 

IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/THE/AN, DF, DS, RED 
PIE=3.14159265359D00 

ALPHA=2 .0D0*DSQORT(S) / (AN+1.0D0-DS* (AN-1.0D0) ) 
BETA=( (DF/DS) * (AN+1.0D0) -DF* (AN-1.0D0) )/2.0D00 
V1=1.00D00-DS/ (AN+1.0D0-DS* (AN-1.0D0) ) 
V2=1.00D00-V1 

ARGU1=ALPHA* (RED**BETA) 

ARGU2=V2*PIE 


NN=1 

CALL DBSKES (V2, ARGU1,NN, BK1) 

CALL DBSIES (V2, ARGU1,NN,BI1) 

T1=BI1/DEXP (ARGU1) +2.0D0O*DSIN (ARGU2) *BK1/PIE/DEXP (ARGU1) 
T6=T1l 

CALL DBSKES (V1, ALPHA, NN, BK2) 
T2=BK2 /DEXP (ALPHA) 

CALL DBSIES (V1, ALPHA, NN, BI3) 

T3=B13/DEXP (ALPHA) 

T4=BK1/DEXP (ARGU1) 

T7=T4 

CALL DBSKES (V2, ALPHA, NN, BK5) 
T5=BK5 /DEXP (ALPHA) 

CALL DBSIES (V2, ALPHA, NN, BI8) 

T8=B1I8/DEXP (ALPHA) +2.0D0*DSIN(ARGU2) *T5/PIE 


CALCULATE THE PRESSURE SOLUTION, AS GIVEN IN EQUATION (7-27) 
PWDL= (T1*T2+T3*T4) / (T5*T6-T7*T8) /S/DSQRT (S) 
PDPL=S*PWDL 


RETURN 
END 
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***x CONSTANT PRESSURE OUTER BOUNDARY CASE Ea ES tS let tak tabs EB 


SUBROUTINE LAP HAS THE FUNCTIONS THAT ARE USED TO 
CALCULATE THE DIMENSIONLESS PRESSURE AND PRESSURE 
DERIVATIVE FUNCTIONS IN TERMS OF THE LAPLACE VARIABLE, 
S, AND THE SYSTEM PARAMETERS. 


SUBROUTINE LAP (S,PWDL, PDPL) 

IMPLICIT REAL*8 (A-H, O-Z) 

COMMON/THE/AN, DF, DS, RED 
PIE=3.14159265359D00 

ALPHA=2 .0D0*DSQRT(S) / (AN+1.0D0-DS* (AN-1.0D0) ) 
BETA=( (DF/DS) * (AN+1.0D0) -DF* (AN-1.0D0) ) /2.0D00 
V1=1 .00D00-DS/ (AN+1.0D0-DS* (AN-1.0D0) ) 
V2=1.00D00-V1 

ARGU1=ALPHA* (RED**BETA) 

ARGU2=V2*PIE 


NN=1 

CALL DBSIES (V1, ARGU1,NN,BI1) 
T1=BI1/DEXP (ARGU1) 

T6=T 1 

CALL DBSKES (V1, ALPHA, NN, BK2) 
T2=BK2 /DEXP (ALPHA) 

CALL DBSIES (V1, ALPHA, NN, BI3) 
T3=B1I3/DEXP (ALPHA) 

CALL DBSKES (V1, ARGU1, NN, BK4) 
T4=BK4/DEXP (ARGU1) 

T7=14 

CALL DBSKES (V2, ALPHA, NN, BK5) 
T5=BK5 /DEXP (ALPHA) 

CALL DBSIES (V2, ALPHA, NN, BI8) 
T8=B18/DEXP (ALPHA) +2 .0D0*DSIN (ARGU2) *T5/PIE 


CALCULATE THE PRESSURE SOLUTION, AS GIVEN IN EQUATION (7-27) 
PWDL= (T1*T2-T3*T4) / (TS*T6+T7*T8) /S/DSQRT (S) 


PDPL=S*PWDL 
RETURN 


END 
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